Coverage for ibllib/pipes/mesoscope_tasks.py: 82%

1"""The mesoscope data extraction pipeline. 

2 

3The mesoscope pipeline currently comprises raw image registration and timestamps extraction. In 

4the future there will be compression (and potential cropping), FOV metadata extraction, and ROI 

5extraction. 

6 

7Pipeline: 

8 1. Register reference images and upload snapshots and notes to Alyx 

9 2. Run ROI cell detection 

10 3. Calculate the pixel and ROI brain locations and register fields of view to Alyx 

11 4. Compress the raw imaging data 

12 5. Extract the imaging times from the main DAQ 

13""" 

14import json (empty)

15import logging (empty)

16import subprocess (empty)

17import shutil (empty)

18from pathlib import Path (empty)

19from itertools import chain (empty)

20from collections import defaultdict, Counter (empty)

21from fnmatch import fnmatch (empty)

22import enum (empty)

23import re (empty)

24import time (empty)

25 

26import numba as nb (empty)

27import numpy as np (empty)

28import pandas as pd (empty)

29import sparse (empty)

30from scipy.io import loadmat (empty)

31from scipy.interpolate import interpn (empty)

32import one.alf.io as alfio (empty)

33from one.alf.spec import is_valid, to_alf (empty)

34from one.alf.files import filename_parts, session_path_parts (empty)

35import one.alf.exceptions as alferr (empty)

36 

37from ibllib.pipes import base_tasks (empty)

38from ibllib.io.extractors import mesoscope (empty)

39from iblatlas.atlas import ALLEN_CCF_LANDMARKS_MLAPDV_UM, MRITorontoAtlas (empty)

40 

41 

42_logger = logging.getLogger(__name__) (empty)

43Provenance = enum.Enum('Provenance', ['ESTIMATE', 'FUNCTIONAL', 'LANDMARK', 'HISTOLOGY']) # py3.11 make StrEnum (empty)

44 

45 

46class MesoscopeRegisterSnapshots(base_tasks.MesoscopeTask, base_tasks.RegisterRawDataTask): (empty)

47 """Upload snapshots as Alyx notes and register the 2P reference image(s).""" 

48 priority = 100 (empty)

49 job_size = 'small' (empty)

50 

51 @property (empty)

52 def signature(self): (empty)

53 signature = { 2 ctx1ol

54 'input_files': [('referenceImage.raw.tif', f'{self.device_collection}/reference', False), 

55 ('referenceImage.stack.tif', f'{self.device_collection}/reference', False), 

56 ('referenceImage.meta.json', f'{self.device_collection}/reference', False)], 

57 'output_files': [('referenceImage.raw.tif', f'{self.device_collection}/reference', False), 

58 ('referenceImage.stack.tif', f'{self.device_collection}/reference', False), 

59 ('referenceImage.meta.json', f'{self.device_collection}/reference', False)] 

60 } 

61 return signature 2 ctx1ol

62 

63 def __init__(self, session_path, **kwargs): (empty)

64 super().__init__(session_path, **kwargs) 4 ctx1polq

65 self.device_collection = self.get_device_collection('mesoscope', 4 ctx1polq

66 kwargs.get('device_collection', 'raw_imaging_data_*')) 

67 

68 def _run(self): (empty)

69 """ 

70 Assert one reference image per collection and rename it. Register snapshots. 

71 

72 Returns 

73 ------- 

74 list of pathlib.Path containing renamed reference image. 

75 """ 

76 # Assert that only one tif file exists per collection 

77 file, collection, _ = self.signature['input_files'][0] 1 ctx1l

78 reference_images = list(self.session_path.rglob(f'{collection}/{file}')) 1 ctx1l

79 assert len(set(x.parent for x in reference_images)) == len(reference_images) 1 ctx1l

80 # Rename the reference images 

81 out_files = super()._run() 1 ctx1l

82 # Register snapshots in base session folder and raw_imaging_data folders 

83 self.register_snapshots(collection=[self.device_collection, '']) 1 ctx1l

84 return out_files 1 ctx1l

85 

86 

87class MesoscopeCompress(base_tasks.MesoscopeTask): (empty)

88 """ Tar compress raw 2p tif files, optionally remove uncompressed data.""" 

89 

90 priority = 90 (empty)

91 job_size = 'large' (empty)

92 _log_level = None (empty)

93 

94 @property (empty)

95 def signature(self): (empty)

96 signature = { 1 ctx1e

97 'input_files': [('*.tif', self.device_collection, True)], 

98 'output_files': [('imaging.frames.tar.bz2', self.device_collection, True)] 

99 } 

100 return signature 1 ctx1e

101 

102 def setUp(self, **kwargs): (empty)

103 """Run at higher log level""" 

104 self._log_level = _logger.level 1 ctx1e

105 _logger.setLevel(logging.DEBUG) 1 ctx1e

106 return super().setUp(**kwargs) 1 ctx1e

107 

108 def tearDown(self): (empty)

109 _logger.setLevel(self._log_level or logging.INFO) 1 ctx1e

110 return super().tearDown() 1 ctx1e

111 

112 def _run(self, remove_uncompressed=False, verify_output=True, clobber=False, **kwargs): (empty)

113 """ 

114 Run tar compression on all tif files in the device collection. 

115 

116 Parameters 

117 ---------- 

118 remove_uncompressed: bool 

119 Whether to remove the original, uncompressed data. Default is False. 

120 verify_output: bool 

121 Whether to check that the compressed tar file can be uncompressed without errors. 

122 Default is True. 

123 

124 Returns 

125 ------- 

126 list of pathlib.Path 

127 Path to compressed tar file. 

128 """ 

129 outfiles = [] # should be one per raw_imaging_data folder 1 ctx1e

130 input_files = defaultdict(list) 1 ctx1e

131 for file, in_dir, _ in self.input_files: 1 ctx1e

132 input_files[self.session_path.joinpath(in_dir)].append(file) 1 ctx1e

133 

134 for in_dir, files in input_files.items(): 1 ctx1e

135 outfile = in_dir / self.output_files[0][0] 1 ctx1e

136 if outfile.exists() and not clobber: 1 ctx1e

137 _logger.info('%s already exists; skipping...', outfile.relative_to(self.session_path)) 

138 continue 

139 # glob for all input patterns 

140 infiles = list(chain(*map(lambda x: in_dir.glob(x), files))) 1 ctx1e

141 if not infiles: 1 ctx1e

142 _logger.info('No image files found in %s', in_dir.relative_to(self.session_path)) 

143 continue 

144 

145 _logger.debug( 1 ctx1e

146 'Input files:\n\t%s', '\n\t'.join(map(Path.as_posix, (x.relative_to(self.session_path) for x in infiles))) 

147 ) 

148 

149 uncompressed_size = sum(x.stat().st_size for x in infiles) 1 ctx1e

150 _logger.info('Compressing %i file(s)', len(infiles)) 1 ctx1e

151 cmd = 'tar -cjvf "{output}" "{input}"'.format( 1 ctx1e

152 output=outfile.relative_to(in_dir), input='" "'.join(str(x.relative_to(in_dir)) for x in infiles)) 

153 _logger.debug(cmd) 1 ctx1e

154 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1 ctx1e

155 info, error = process.communicate() # b'2023-02-17_2_test_2P_00001_00001.tif\n' 1 ctx1e

156 _logger.debug(info.decode()) 1 ctx1e

157 assert process.returncode == 0, f'compression failed: {error.decode()}' 1 ctx1e

158 

159 # Check the output 

160 assert outfile.exists(), 'output file missing' 1 ctx1e

161 outfiles.append(outfile) 1 ctx1e

162 compressed_size = outfile.stat().st_size 1 ctx1e

163 min_size = kwargs.pop('verify_min_size', 1024) 1 ctx1e

164 assert compressed_size > int(min_size), f'Compressed file < {min_size / 1024:.0f}KB' 1 ctx1e

165 _logger.info('Compression ratio = %.3f, saving %.2f pct (%.2f MB)', 1 ctx1e

166 uncompressed_size / compressed_size, 

167 round((1 - (compressed_size / uncompressed_size)) * 10000) / 100, 

168 (uncompressed_size - compressed_size) / 1024 / 1024) 

169 

170 if verify_output: 1 ctx1e

171 # Test bzip 

172 cmd = f'bzip2 -tv {outfile.relative_to(in_dir)}' 1 ctx1e

173 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1 ctx1e

174 info, error = process.communicate() 1 ctx1e

175 _logger.debug(info.decode()) 1 ctx1e

176 assert process.returncode == 0, f'bzip compression test failed: {error}' 1 ctx1e

177 # Check tar 

178 cmd = f'bunzip2 -dc {outfile.relative_to(in_dir)} | tar -tvf -' 1 ctx1e

179 process = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=in_dir) 1 ctx1e

180 info, error = process.communicate() 1 ctx1e

181 _logger.debug(info.decode()) 1 ctx1e

182 assert process.returncode == 0, 'tarball decompression test failed' 1 ctx1e

183 compressed_files = set(x.split()[-1] for x in filter(None, info.decode().split('\n'))) 1 ctx1e

184 assert compressed_files == set(x.name for x in infiles) 1 ctx1e

185 

186 if remove_uncompressed: 1 ctx1e

187 _logger.info(f'Removing input files for {in_dir.relative_to(self.session_path)}') 1 ctx1e

188 for file in infiles: 1 ctx1e

189 file.unlink() 1 ctx1e

190 

191 return outfiles 1 ctx1e

192 

193 

194class MesoscopePreprocess(base_tasks.MesoscopeTask): (empty)

195 """Run suite2p preprocessing on tif files""" 

196 

197 priority = 80 (empty)

198 cpu = -1 (empty)

199 job_size = 'large' (empty)

200 

201 @property (empty)

202 def signature(self): (empty)

203 # The number of in and outputs will be dependent on the number of input raw imaging folders and output FOVs 

204 signature = { 1 ctx1f

205 'input_files': [('_ibl_rawImagingData.meta.json', self.device_collection, True), 

206 ('*.tif', self.device_collection, True), 

207 ('exptQC.mat', self.device_collection, False)], 

208 'output_files': [('mpci.ROIActivityF.npy', 'alf/FOV*', True), 

209 ('mpci.ROINeuropilActivityF.npy', 'alf/FOV*', True), 

210 ('mpci.ROIActivityDeconvolved.npy', 'alf/FOV*', True), 

211 ('mpci.badFrames.npy', 'alf/FOV*', True), 

212 ('mpci.mpciFrameQC.npy', 'alf/FOV*', True), 

213 ('mpciFrameQC.names.tsv', 'alf/FOV*', True), 

214 ('mpciMeanImage.images.npy', 'alf/FOV*', True), 

215 ('mpciROIs.stackPos.npy', 'alf/FOV*', True), 

216 ('mpciROIs.mpciROITypes.npy', 'alf/FOV*', True), 

217 ('mpciROIs.cellClassifier.npy', 'alf/FOV*', True), 

218 ('mpciROITypes.names.tsv', 'alf/FOV*', True), 

219 ('mpciROIs.masks.npy', 'alf/FOV*', True), 

220 ('mpciROIs.neuropilMasks.npy', 'alf/FOV*', True), 

221 ('_suite2p_ROIData.raw.zip', self.device_collection, False)] 

222 } 

223 return signature 1 ctx1f

224 

225 @staticmethod (empty)

226 def _masks2sparse(stat, ops): (empty)

227 """ 

228 Extract 3D sparse mask arrays from suit2p output. 

229 

230 Parameters 

231 ---------- 

232 stat : numpy.array 

233 The loaded stat.npy file. A structured array with fields ('lam', 'ypix', 'xpix', 'neuropil_mask'). 

234 ops : numpy.array 

235 The loaded ops.npy file. A structured array with fields ('Ly', 'Lx'). 

236 

237 Returns 

238 ------- 

239 sparse.GCXS 

240 A pydata sparse array of type float32, representing the ROI masks. 

241 sparse.GCXS 

242 A pydata sparse array of type float32, representing the neuropil ROI masks. 

243 

244 Notes 

245 ----- 

246 Save using sparse.save_npz. 

247 """ 

248 shape = (stat.shape[0], ops['Ly'], ops['Lx']) 2 ctx1dc

249 npx = np.prod(shape[1:]) # Number of pixels per time point 2 ctx1dc

250 coords = [[], [], []] 2 ctx1dc

251 data = [] 2 ctx1dc

252 pil_coords = [] 2 ctx1dc

253 for i, s in enumerate(stat): 2 ctx1dc

254 coords[0].append(np.full(s['ypix'].shape, i)) 2 ctx1dc

255 coords[1].append(s['ypix']) 2 ctx1dc

256 coords[2].append(s['xpix']) 2 ctx1dc

257 data.append(s['lam']) 2 ctx1dc

258 pil_coords.append(s['neuropil_mask'] + i * npx) 2 ctx1dc

259 roi_mask_sp = sparse.COO(list(map(np.concatenate, coords)), np.concatenate(data), shape=shape) 2 ctx1dc

260 pil_mask_sp = sparse.COO(np.unravel_index(np.concatenate(pil_coords), shape), True, shape=shape) 2 ctx1dc

261 return sparse.GCXS.from_coo(roi_mask_sp), sparse.GCXS.from_coo(pil_mask_sp) 2 ctx1dc

262 

263 def _rename_outputs(self, suite2p_dir, frameQC_names, frameQC, rename_dict=None): (empty)

264 """ 

265 Convert suite2p output files to ALF datasets. 

266 

267 Parameters 

268 ---------- 

269 suite2p_dir : pathlib.Path 

270 rename_dict : dict or None 

271 The suite2p output filenames and the corresponding ALF name. NB: These files are saved 

272 after transposition. Default is None, i.e. using the default mapping hardcoded in the function below. 

273 

274 Returns 

275 ------- 

276 list of pathlib.Path 

277 All paths found in FOV folders. 

278 """ 

279 if rename_dict is None: 2 ctx1dc

280 rename_dict = { 2 ctx1dc

281 'F.npy': 'mpci.ROIActivityF.npy', 

282 'spks.npy': 'mpci.ROIActivityDeconvolved.npy', 

283 'Fneu.npy': 'mpci.ROINeuropilActivityF.npy' 

284 } 

285 # Rename the outputs, first the subdirectories 

286 for plane_dir in suite2p_dir.iterdir(): 2 ctx1dc

287 # ignore the combined dir 

288 if plane_dir.name != 'combined': 2 ctx1dc

289 n = int(plane_dir.name.split('plane')[1]) 2 ctx1dc

290 fov_dir = plane_dir.parent.joinpath(f'FOV_{n:02}') 2 ctx1dc

291 if fov_dir.exists(): 2 ctx1dc

292 shutil.rmtree(str(fov_dir), ignore_errors=False, onerror=None) 

293 plane_dir.rename(fov_dir) 2 ctx1dc

294 # Now rename the content of the new directories and move them out of suite2p 

295 for fov_dir in suite2p_dir.iterdir(): 2 ctx1dc

296 # Compress suite2p output files 

297 target = suite2p_dir.parent.joinpath(fov_dir.name) 2 ctx1dc

298 target.mkdir(exist_ok=True) 2 ctx1dc

299 # Move bin file out of the way first 

300 if fov_dir.joinpath('data.bin').exists(): 2 ctx1dc

301 dst = self.session_path.joinpath('raw_bin_files', fov_dir.name, 'data.bin') 

302 dst.parent.mkdir(parents=True, exist_ok=True) 

303 _logger.debug('Moving bin file to %s', dst.relative_to(self.session_path)) 

304 fov_dir.joinpath('data.bin').replace(dst) 

305 # Set logger to warning for the moment to not clutter the logs 

306 prev_level = _logger.level 2 ctx1dc

307 _logger.setLevel(logging.WARNING) 2 ctx1dc

308 shutil.make_archive(str(target / '_suite2p_ROIData.raw'), 'zip', fov_dir, logger=_logger) 2 ctx1dc

309 _logger.setLevel(prev_level) 2 ctx1dc

310 if fov_dir != 'combined': 2 ctx1dc

311 # save frameQC in each dir (for now, maybe there will be fov specific frame QC eventually) 

312 if frameQC is not None and len(frameQC) > 0: 2 ctx1dc

313 np.save(fov_dir.joinpath('mpci.mpciFrameQC.npy'), frameQC) 1 ctx1c

314 frameQC_names.to_csv(fov_dir.joinpath('mpciFrameQC.names.tsv'), sep='\t', index=False) 1 ctx1c

315 

316 # extract some other data from suite2p outputs 

317 ops = np.load(fov_dir.joinpath('ops.npy'), allow_pickle=True).item() 2 ctx1dc

318 stat = np.load(fov_dir.joinpath('stat.npy'), allow_pickle=True) 2 ctx1dc

319 iscell = np.load(fov_dir.joinpath('iscell.npy')) 2 ctx1dc

320 

321 # Save suite2p ROI activity outputs in transposed from (n_frames, n_ROI) 

322 for k, v in rename_dict.items(): 2 ctx1dc

323 np.save(fov_dir.joinpath(v), np.load(fov_dir.joinpath(k)).T) 2 ctx1dc

324 # fov_dir.joinpath(k).unlink() # Keep original files for suite2P GUI 

325 np.save(fov_dir.joinpath('mpci.badFrames.npy'), np.asarray(ops['badframes'], dtype=bool)) 2 ctx1dc

326 np.save(fov_dir.joinpath('mpciMeanImage.images.npy'), np.asarray(ops['meanImg'], dtype=float)) 2 ctx1dc

327 np.save(fov_dir.joinpath('mpciROIs.stackPos.npy'), np.asarray([(*s['med'], 0) for s in stat], dtype=int)) 2 ctx1dc

328 np.save(fov_dir.joinpath('mpciROIs.cellClassifier.npy'), np.asarray(iscell[:, 1], dtype=float)) 2 ctx1dc

329 np.save(fov_dir.joinpath('mpciROIs.mpciROITypes.npy'), np.asarray(iscell[:, 0], dtype=np.int16)) 2 ctx1dc

330 pd.DataFrame([(0, 'no cell'), (1, 'cell')], columns=['roi_values', 'roi_labels'] 2 ctx1dc

331 ).to_csv(fov_dir.joinpath('mpciROITypes.names.tsv'), sep='\t', index=False) 

332 # ROI and neuropil masks 

333 roi_mask, pil_mask = self._masks2sparse(stat, ops) 2 ctx1dc

334 with open(fov_dir.joinpath('mpciROIs.masks.sparse_npz'), 'wb') as fp: 2 ctx1dc

335 sparse.save_npz(fp, roi_mask) 2 ctx1dc

336 with open(fov_dir.joinpath('mpciROIs.neuropilMasks.sparse_npz'), 'wb') as fp: 2 ctx1dc

337 sparse.save_npz(fp, pil_mask) 2 ctx1dc

338 # move folders out of suite2p dir 

339 # We overwrite existing files 

340 for file in filter(lambda x: is_valid(x.name), fov_dir.iterdir()): 2 ctx1dc

341 target_file = target.joinpath(file.name) 2 ctx1dc

342 if target_file.exists(): 2 ctx1dc

343 target_file.unlink() 

344 file.rename(target_file) 2 ctx1dc

345 shutil.rmtree(str(suite2p_dir), ignore_errors=False, onerror=None) 2 ctx1dc

346 # Collect all files in those directories 

347 return list(suite2p_dir.parent.rglob('FOV*/*')) 2 ctx1dc

348 

349 @staticmethod (empty)

350 def _check_meta_data(meta_data_all: list) -> dict: (empty)

351 """ 

352 Check that the meta data is consistent across all raw imaging folders. 

353 

354 Parameters 

355 ---------- 

356 meta_data_all: list of dicts 

357 List of metadata dictionaries to be checked for consistency. 

358 

359 Returns 

360 ------- 

361 dict 

362 Single, consolidated dictionary containing metadata. 

363 """ 

364 # Ignore the things we don't expect to match 

365 ignore = ('acquisitionStartTime', 'nFrames') 

366 ignore_sub = {'rawScanImageMeta': ('ImageDescription', 'Software')} 

367 

368 def equal_dicts(a, b, skip=None): 

369 ka = set(a).difference(skip or ()) 

370 kb = set(b).difference(skip or ()) 

371 return ka == kb and all(a[key] == b[key] for key in ka) 

372 

373 # Compare each dict with the first one in the list 

374 for i, meta in enumerate(meta_data_all[1:]): 

375 if meta != meta_data_all[0]: # compare entire object first 

376 for k, v in meta_data_all[0].items(): # check key by key 

377 if not (equal_dicts(v, meta[k], ignore_sub[k]) # compare sub-dicts... 

378 if k in ignore_sub else # ... if we have keys to ignore in test 

379 not (k in ignore or v == meta[k])): 

380 _logger.warning(f'Mismatch in meta data between raw_imaging_data folders for key {k}. ' 

381 f'Using meta_data from first folder!') 

382 return meta_data_all[0] 

383 

384 @staticmethod (empty)

385 def _consolidate_exptQC(exptQC): (empty)

386 """ 

387 Consolidate exptQC.mat files into a single file. 

388 

389 Parameters 

390 ---------- 

391 exptQC : list of pandas.DataFrame 

392 The loaded 'exptQC.mat' files as squeezed and simplified data frames, with columns 

393 {'frameQC_frames', 'frameQC_names'}. 

394 

395 Returns 

396 ------- 

397 numpy.array 

398 An array of uint8 where 0 indicates good frames, and other values correspond to 

399 experimenter-defined issues (in 'qc_values' column of output data frame). 

400 pandas.DataFrame 

401 A data frame with columns {'qc_values', 'qc_labels'}, the former an unsigned int 

402 corresponding to a QC code; the latter a human-readable QC explanation. 

403 numpy.array 

404 An array of frame indices where QC code != 0. 

405 """ 

406 

407 # Merge and make sure same indexes have same names across all files 

408 frameQC_names_list = [e['frameQC_names'] for e in exptQC] 

409 frameQC_names_list = [{f: 0} if isinstance(f, str) else {f[i]: i for i in range(len(f))} 

410 for f in frameQC_names_list] 

411 frameQC_names = {k: v for d in frameQC_names_list for k, v in d.items()} 

412 for d in frameQC_names_list: 

413 for k, v in d.items(): 

414 if frameQC_names[k] != v: 

415 raise IOError(f'exptQC.mat files have different values for name "{k}"') 

416 

417 frameQC_names = pd.DataFrame(sorted([(v, k) for k, v in frameQC_names.items()]), 

418 columns=['qc_values', 'qc_labels']) 

419 

420 # Concatenate frames 

421 frameQC = np.concatenate([e['frameQC_frames'] for e in exptQC], axis=0) 

422 

423 # Transform to bad_frames as expected by suite2p 

424 bad_frames = np.where(frameQC != 0)[0] 

425 

426 return frameQC, frameQC_names, bad_frames 

427 

428 def get_default_tau(self): (empty)

429 """ 

430 Determine the tau (fluorescence decay) from the subject's genotype. 

431 

432 Returns 

433 ------- 

434 float 

435 The tau value to use. 

436 

437 See Also 

438 -------- 

439 https://suite2p.readthedocs.io/en/latest/settings.html 

440 """ 

441 # These settings are from the suite2P documentation 

442 TAU_MAP = {'G6s': 1.5, 'G6m': 1., 'G6f': .7, 'default': 1.5} 1 ctx1m

443 _, subject, *_ = session_path_parts(self.session_path) 1 ctx1m

444 genotype = self.one.alyx.rest('subjects', 'read', id=subject)['genotype'] 1 ctx1m

445 match = next(filter(None, (re.match(r'.+-(G\d[fms])$', g['allele']) for g in genotype)), None) 1 ctx1m

446 key = match.groups()[0] if match else 'default' 1 ctx1m

447 return TAU_MAP.get(key, TAU_MAP['default']) 1 ctx1m

448 

449 def _create_db(self, meta): (empty)

450 """ 

451 Create the ops dictionary for suite2p based on metadata. 

452 

453 Parameters 

454 ---------- 

455 meta: dict 

456 Imaging metadata. 

457 

458 Returns 

459 ------- 

460 dict 

461 Inputs to suite2p run that deviate from default parameters. 

462 """ 

463 

464 # Currently only supporting single plane, assert that this is the case 

465 # FIXME This checks for zstacks but not dual plane mode 

466 if not isinstance(meta['scanImageParams']['hStackManager']['zs'], int): 1 ctx1f

467 raise NotImplementedError('Multi-plane imaging not yet supported, data seems to be multi-plane') 

468 

469 # Computing dx and dy 

470 cXY = np.array([fov['topLeftDeg'] for fov in meta['FOV']]) 1 ctx1f

471 cXY -= np.min(cXY, axis=0) 1 ctx1f

472 nXnYnZ = np.array([fov['nXnYnZ'] for fov in meta['FOV']]) 1 ctx1f

473 sW = np.sqrt(np.sum((np.array([fov['topRightDeg'] for fov in meta['FOV']]) - np.array( 1 ctx1f

474 [fov['topLeftDeg'] for fov in meta['FOV']])) ** 2, axis=1)) 

475 sH = np.sqrt(np.sum((np.array([fov['bottomLeftDeg'] for fov in meta['FOV']]) - np.array( 1 ctx1f

476 [fov['topLeftDeg'] for fov in meta['FOV']])) ** 2, axis=1)) 

477 pixSizeX = nXnYnZ[:, 0] / sW 1 ctx1f

478 pixSizeY = nXnYnZ[:, 1] / sH 1 ctx1f

479 dx = np.round(cXY[:, 0] * pixSizeX).astype(dtype=np.int32) 1 ctx1f

480 dy = np.round(cXY[:, 1] * pixSizeY).astype(dtype=np.int32) 1 ctx1f

481 nchannels = len(meta['channelSaved']) if isinstance(meta['channelSaved'], list) else 1 1 ctx1f

482 

483 db = { 1 ctx1f

484 'data_path': sorted(map(str, self.session_path.glob(f'{self.device_collection}'))), 

485 'save_path0': str(self.session_path.joinpath('alf')), 

486 'fast_disk': '', # TODO 

487 'look_one_level_down': False, # don't look in the children folders as that is where the reference data is 

488 'num_workers': self.cpu, # this selects number of cores to parallelize over for the registration step 

489 'num_workers_roi': -1, # for parallelization over FOVs during cell detection, for now don't 

490 'keep_movie_raw': False, 

491 'delete_bin': False, # TODO: delete this on the long run 

492 'batch_size': 500, # SP reduced this from 1000 

493 'nimg_init': 400, 

494 'combined': False, 

495 'nonrigid': True, 

496 'maxregshift': 0.05, # default = 1 

497 'denoise': 1, # whether binned movie should be denoised before cell detection 

498 'block_size': [128, 128], 

499 'save_mat': True, # save the data to Fall.mat 

500 'move_bin': True, # move the binary file to save_path 

501 'scalefactor': 1, # scale manually in x to account for overlap between adjacent ribbons UCL mesoscope 

502 'mesoscan': True, 

503 'nplanes': 1, 

504 'nrois': len(meta['FOV']), 

505 'nchannels': nchannels, 

506 'fs': meta['scanImageParams']['hRoiManager']['scanVolumeRate'], 

507 'lines': [list(np.asarray(fov['lineIdx']) - 1) for fov in meta['FOV']], # subtracting 1 to make 0-based 

508 'tau': self.get_default_tau(), # deduce the GCamp used from Alyx mouse line (defaults to 1.5; that of GCaMP6s) 

509 'functional_chan': 1, # for now, eventually find(ismember(meta.channelSaved == meta.channelID.green)) 

510 'align_by_chan': 1, # for now, eventually find(ismember(meta.channelSaved == meta.channelID.red)) 

511 'dx': dx, 

512 'dy': dy 

513 } 

514 

515 return db 1 ctx1f

516 

517 def _run(self, run_suite2p=True, rename_files=True, use_badframes=True, **kwargs): (empty)

518 """ 

519 Process inputs, run suite2p and make outputs alf compatible. 

520 

521 Parameters 

522 ---------- 

523 run_suite2p: bool 

524 Whether to run suite2p, default is True. 

525 rename_files: bool 

526 Whether to rename and reorganize the suite2p outputs to be alf compatible. Defaults is True. 

527 use_badframes: bool 

528 Whether to exclude bad frames indicated by the experimenter in exptQC.mat. Default is currently False 

529 due to bug in suite2p. Change this in the future. 

530 

531 Returns 

532 ------- 

533 list of pathlib.Path 

534 All files created by the task. 

535 """ 

536 import suite2p 1 ctx1f

537 

538 """ Metadata and parameters """ 1 ctx1f

539 # Load metadata and make sure all metadata is consistent across FOVs 

540 meta_files = sorted(self.session_path.glob(f'{self.device_collection}/*rawImagingData.meta.*')) 1 ctx1f

541 collections = set(f.parts[-2] for f in meta_files) 1 ctx1f

542 # Check there is exactly 1 meta file per collection 

543 assert len(meta_files) == len(list(self.session_path.glob(self.device_collection))) == len(collections) 1 ctx1f

544 rawImagingData = [mesoscope.patch_imaging_meta(alfio.load_file_content(filepath)) for filepath in meta_files] 1 ctx1f

545 if len(rawImagingData) > 1: 1 ctx1f

546 meta = self._check_meta_data(rawImagingData) 

547 else: 

548 meta = rawImagingData[0] 1 ctx1f

549 # Get default ops 

550 ops = suite2p.default_ops() 1 ctx1f

551 # Create db which overwrites ops when passed to suite2p, with information from meta data and hardcoded 

552 db = self._create_db(meta) 1 ctx1f

553 # Anything can be overwritten by keyword arguments passed to the tasks run() method 

554 for k, v in kwargs.items(): 1 ctx1f

555 if k in ops.keys() or k in db.keys(): 1 ctx1f

556 # db overwrites ops when passed to run_s2p, so we only need to update / add it here 

557 db[k] = v 1 ctx1f

558 # Update the task kwargs attribute as it will be stored in the arguments json field in alyx 

559 self.kwargs = {**self.kwargs, **db} 1 ctx1f

560 

561 """ Bad frames """ 1 ctx1f

562 qc_paths = (self.session_path.joinpath(f[1], 'exptQC.mat') 1 ctx1f

563 for f in self.input_files if f[0] == 'exptQC.mat') 

564 qc_paths = map(str, filter(Path.exists, qc_paths)) 1 ctx1f

565 exptQC = [loadmat(p, squeeze_me=True, simplify_cells=True) for p in qc_paths] 1 ctx1f

566 if len(exptQC) > 0: 1 ctx1f

567 frameQC, frameQC_names, bad_frames = self._consolidate_exptQC(exptQC) 

568 else: 

569 _logger.warning('No frame QC (exptQC.mat) files found.') 1 ctx1f

570 frameQC, bad_frames = np.array([], dtype='u1'), np.array([], dtype='i8') 1 ctx1f

571 frameQC_names = pd.DataFrame(columns=['qc_values', 'qc_labels']) 1 ctx1f

572 # If applicable, save as bad_frames.npy in first raw_imaging_folder for suite2p 

573 if len(bad_frames) > 0 and use_badframes is True: 1 ctx1f

574 np.save(Path(db['data_path'][0]).joinpath('bad_frames.npy'), bad_frames) 

575 

576 """ Suite2p """ 1 ctx1f

577 # Create alf it is doesn't exist 

578 self.session_path.joinpath('alf').mkdir(exist_ok=True) 1 ctx1f

579 # Remove existing suite2p dir if it exists 

580 suite2p_dir = Path(db['save_path0']).joinpath('suite2p') 1 ctx1f

581 if suite2p_dir.exists(): 1 ctx1f

582 shutil.rmtree(str(suite2p_dir), ignore_errors=True, onerror=None) 

583 # Run suite2p 

584 if run_suite2p: 1 ctx1f

585 _ = suite2p.run_s2p(ops=ops, db=db) 

586 

587 """ Outputs """ 1 ctx1f

588 # Save and rename other outputs 

589 if rename_files: 1 ctx1f

590 out_files = self._rename_outputs(suite2p_dir, frameQC_names, frameQC) 

591 else: 

592 out_files = list(Path(db['save_path0']).joinpath('suite2p').rglob('*')) 1 ctx1f

593 # Only return output file that are in the signature (for registration) 

594 out_files = [f for f in out_files if f.name in [f[0] for f in self.output_files]] 1 ctx1f

595 return out_files 1 ctx1f

596 

597 

598class MesoscopeSync(base_tasks.MesoscopeTask): (empty)

599 """Extract the frame times from the main DAQ.""" 

600 

601 priority = 40 (empty)

602 job_size = 'small' (empty)

603 

604 @property (empty)

605 def signature(self): (empty)

606 signature = { 2 ctx1ij

607 'input_files': [(f'_{self.sync_namespace}_DAQdata.raw.npy', self.sync_collection, True), 

608 (f'_{self.sync_namespace}_DAQdata.timestamps.npy', self.sync_collection, True), 

609 (f'_{self.sync_namespace}_DAQdata.meta.json', self.sync_collection, True), 

610 ('_ibl_rawImagingData.meta.json', self.device_collection, True), 

611 ('rawImagingData.times_scanImage.npy', self.device_collection, True), 

612 (f'_{self.sync_namespace}_softwareEvents.log.htsv', self.sync_collection, False), ], 

613 'output_files': [('mpci.times.npy', 'alf/mesoscope/FOV*', True), 

614 ('mpciStack.timeshift.npy', 'alf/mesoscope/FOV*', True), ] 

615 } 

616 return signature 2 ctx1ij

617 

618 def _run(self): (empty)

619 """ 

620 Extract the imaging times for all FOVs. 

621 

622 Returns 

623 ------- 

624 list of pathlib.Path 

625 Files containing frame timestamps for individual FOVs and time offsets for each line scan. 

626 

627 """ 

628 # TODO function to determine nFOVs 

629 try: 2 ctx1ij

630 alf_path = self.session_path / self.sync_collection 2 ctx1ij

631 events = alfio.load_object(alf_path, 'softwareEvents').get('log') 2 ctx1ij

632 except alferr.ALFObjectNotFound: 1 ctx1i

633 events = None 1 ctx1i

634 if events is None or events.empty: 2 ctx1ij

635 _logger.debug('No software events found for session %s', self.session_path) 1 ctx1i

636 collections = set(collection for _, collection, _ in self.input_files 2 ctx1ij

637 if fnmatch(collection, self.device_collection)) 

638 # Load first meta data file to determine the number of FOVs 

639 # Changing FOV between imaging bouts is not supported currently! 

640 self.rawImagingData = alfio.load_object(self.session_path / next(iter(collections)), 'rawImagingData') 2 ctx1ij

641 self.rawImagingData['meta'] = mesoscope.patch_imaging_meta(self.rawImagingData['meta']) 2 ctx1ij

642 n_FOVs = len(self.rawImagingData['meta']['FOV']) 2 ctx1ij

643 sync, chmap = self.load_sync() # Extract sync data from raw DAQ data 2 ctx1ij

644 mesosync = mesoscope.MesoscopeSyncTimeline(self.session_path, n_FOVs) 2 ctx1ij

645 _, out_files = mesosync.extract( 2 ctx1ij

646 save=True, sync=sync, chmap=chmap, device_collection=collections, events=events) 

647 return out_files 2 ctx1ij

648 

649 

650class MesoscopeFOV(base_tasks.MesoscopeTask): (empty)

651 """Create FOV and FOV location objects in Alyx from metadata.""" 

652 

653 priority = 40 (empty)

654 job_size = 'small' (empty)