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

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)

18import uuid (empty)

19from pathlib import Path (empty)

20from itertools import chain (empty)

21from collections import defaultdict, Counter (empty)

22from fnmatch import fnmatch (empty)

23import enum (empty)

24import re (empty)

25import time (empty)

26 

27import numba as nb (empty)

28import numpy as np (empty)

29import pandas as pd (empty)

30import sparse (empty)

31from scipy.io import loadmat (empty)

32from scipy.interpolate import interpn (empty)

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

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

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

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

37 

38from ibllib.pipes import base_tasks (empty)

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

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

41 

42 

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

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

45 

46 

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

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

49 priority = 100 (empty)

50 job_size = 'small' (empty)

51 

52 @property (empty)

53 def signature(self): (empty)

54 signature = { 2 ctx1pm

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

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

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

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

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

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

61 } 

62 return signature 2 ctx1pm

63 

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

65 super().__init__(session_path, **kwargs) 3 ctx1qpm

66 self.device_collection = self.get_device_collection('mesoscope', 3 ctx1qpm

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

68 

69 def _run(self): (empty)

70 """ 

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

72 

73 Returns 

74 ------- 

75 list of pathlib.Path containing renamed reference image. 

76 """ 

77 # Assert that only one tif file exists per collection 

78 file, collection, _ = self.signature['input_files'][0] 1 ctx1m

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

80 assert len(set(x.parent for x in reference_images)) == len(reference_images) 1 ctx1m

81 # Rename the reference images 

82 out_files = super()._run() 1 ctx1m

83 # Register snapshots in base session folder and raw_imaging_data folders 

84 self.register_snapshots(collection=[self.device_collection, '']) 1 ctx1m

85 return out_files 1 ctx1m

86 

87 

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

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

90 

91 priority = 90 (empty)

92 job_size = 'large' (empty)

93 _log_level = None (empty)

94 

95 @property (empty)

96 def signature(self): (empty)

97 signature = { 1 ctx1f

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

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

100 } 

101 return signature 1 ctx1f

102 

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

104 """Run at higher log level""" 

105 self._log_level = _logger.level 1 ctx1f

106 _logger.setLevel(logging.DEBUG) 1 ctx1f

107 return super().setUp(**kwargs) 1 ctx1f

108 

109 def tearDown(self): (empty)

110 _logger.setLevel(self._log_level or logging.INFO) 1 ctx1f

111 return super().tearDown() 1 ctx1f

112 

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

114 """ 

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

116 

117 Parameters 

118 ---------- 

119 remove_uncompressed: bool 

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

121 verify_output: bool 

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

123 Default is True. 

124 

125 Returns 

126 ------- 

127 list of pathlib.Path 

128 Path to compressed tar file. 

129 """ 

130 outfiles = [] # should be one per raw_imaging_data folder 1 ctx1f

131 input_files = defaultdict(list) 1 ctx1f

132 for file, in_dir, _ in self.input_files: 1 ctx1f

133 input_files[self.session_path.joinpath(in_dir)].append(file) 1 ctx1f

134 

135 for in_dir, files in input_files.items(): 1 ctx1f

136 outfile = in_dir / self.output_files[0][0] 1 ctx1f

137 if outfile.exists() and not clobber: 1 ctx1f

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

139 continue 

140 # glob for all input patterns 

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

142 if not infiles: 1 ctx1f

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

144 continue 

145 

146 _logger.debug( 1 ctx1f

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

148 ) 

149 

150 uncompressed_size = sum(x.stat().st_size for x in infiles) 1 ctx1f

151 _logger.info('Compressing %i file(s)', len(infiles)) 1 ctx1f

152 cmd = 'tar -cjvf "{output}" "{input}"'.format( 1 ctx1f

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

154 _logger.debug(cmd) 1 ctx1f

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

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

157 _logger.debug(info.decode()) 1 ctx1f

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

159 

160 # Check the output 

161 assert outfile.exists(), 'output file missing' 1 ctx1f

162 outfiles.append(outfile) 1 ctx1f

163 compressed_size = outfile.stat().st_size 1 ctx1f

164 min_size = kwargs.pop('verify_min_size', 1024) 1 ctx1f

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

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

167 uncompressed_size / compressed_size, 

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

169 (uncompressed_size - compressed_size) / 1024 / 1024) 

170 

171 if verify_output: 1 ctx1f

172 # Test bzip 

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

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

175 info, error = process.communicate() 1 ctx1f

176 _logger.debug(info.decode()) 1 ctx1f

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

178 # Check tar 

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

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

181 info, error = process.communicate() 1 ctx1f

182 _logger.debug(info.decode()) 1 ctx1f

183 assert process.returncode == 0, 'tarball decompression test failed' 1 ctx1f

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

185 assert compressed_files == set(x.name for x in infiles) 1 ctx1f

186 

187 if remove_uncompressed: 1 ctx1f

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

189 for file in infiles: 1 ctx1f

190 file.unlink() 1 ctx1f

191 

192 return outfiles 1 ctx1f

193 

194 

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

196 """Run suite2p preprocessing on tif files.""" 

197 

198 priority = 80 (empty)

199 cpu = -1 (empty)

200 job_size = 'large' (empty)

201 env = 'suite2p' (empty)

202 

203 @property (empty)

204 def signature(self): (empty)

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

206 signature = { 1 ctx1e

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

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

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

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

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

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

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

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

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

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

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

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

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

220 ('mpciROIs.uuids.csv', 'alf/FOV*', True), 

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

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

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

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

225 } 

226 return signature 1 ctx1e

227 

228 @staticmethod (empty)

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

230 """ 

231 Extract 3D sparse mask arrays from suit2p output. 

232 

233 Parameters 

234 ---------- 

235 stat : numpy.array 

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

237 ops : numpy.array 

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

239 

240 Returns 

241 ------- 

242 sparse.GCXS 

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

244 sparse.GCXS 

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

246 

247 Notes 

248 ----- 

249 Save using sparse.save_npz. 

250 """ 

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

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

253 coords = [[], [], []] 2 ctx1dc

254 data = [] 2 ctx1dc

255 pil_coords = [] 2 ctx1dc

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

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

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

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

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

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

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

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

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

265 

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

267 """ 

268 Convert suite2p output files to ALF datasets. 

269 

270 Parameters 

271 ---------- 

272 suite2p_dir : pathlib.Path 

273 rename_dict : dict or None 

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

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

276 

277 Returns 

278 ------- 

279 list of pathlib.Path 

280 All paths found in FOV folders. 

281 """ 

282 if rename_dict is None: 2 ctx1dc

283 rename_dict = { 2 ctx1dc

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

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

286 'Fneu.npy': 'mpci.ROINeuropilActivityF.npy' 

287 } 

288 # Rename the outputs, first the subdirectories 

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

290 # ignore the combined dir 

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

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

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

294 if fov_dir.exists(): 2 ctx1dc

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

296 plane_dir.rename(fov_dir) 2 ctx1dc

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

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

299 # Compress suite2p output files 

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

301 target.mkdir(exist_ok=True) 2 ctx1dc

302 # Move bin file out of the way first 

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

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

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

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

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

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

309 prev_level = _logger.level 2 ctx1dc

310 _logger.setLevel(logging.WARNING) 2 ctx1dc

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

312 _logger.setLevel(prev_level) 2 ctx1dc

313 if fov_dir != 'combined': 2 ctx1dc

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

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

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

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

318 

319 # extract some other data from suite2p outputs 

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

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

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

323 

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

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

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

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

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

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

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

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

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

333 # clusters uuids 

334 uuid_list = ['uuids'] + list(map(str, [uuid.uuid4() for _ in range(len(iscell))])) 2 ctx1dc

335 with open(fov_dir.joinpath('mpciROIs.uuids.csv'), 'w+') as fid: 2 ctx1dc

336 fid.write('\n'.join(uuid_list)) 2 ctx1dc

337 

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

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

340 # ROI and neuropil masks 

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

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

343 sparse.save_npz(fp, roi_mask) 2 ctx1dc

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

345 sparse.save_npz(fp, pil_mask) 2 ctx1dc

346 # move folders out of suite2p dir 

347 # We overwrite existing files 

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

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

350 if target_file.exists(): 2 ctx1dc

351 target_file.unlink() 

352 file.rename(target_file) 2 ctx1dc

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

354 # Collect all files in those directories 

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

356 

357 @staticmethod (empty)

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

359 """ 

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

361 

362 Parameters 

363 ---------- 

364 meta_data_all: list of dicts 

365 List of metadata dictionaries to be checked for consistency. 

366 

367 Returns 

368 ------- 

369 dict 

370 Single, consolidated dictionary containing metadata. 

371 """ 

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

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

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

375 

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

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

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

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

380 

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

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

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

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

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

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

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

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

389 f'Using meta_data from first folder!') 

390 return meta_data_all[0] 

391 

392 @staticmethod (empty)

393 def _consolidate_exptQC(exptQC): (empty)

394 """ 

395 Consolidate exptQC.mat files into a single file. 

396 

397 Parameters 

398 ---------- 

399 exptQC : list of pandas.DataFrame 

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

401 {'frameQC_frames', 'frameQC_names'}. 

402 

403 Returns 

404 ------- 

405 numpy.array 

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

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

408 pandas.DataFrame 

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

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

411 numpy.array 

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

413 """ 

414 

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

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

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

418 for f in frameQC_names_list] 

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

420 for d in frameQC_names_list: 

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

422 if frameQC_names[k] != v: 

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

424 

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

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

427 

428 # Concatenate frames 

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

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

431 return frameQC, frameQC_names, bad_frames 

432 

433 def get_default_tau(self): (empty)

434 """ 

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

436 

437 Returns 

438 ------- 

439 float 

440 The tau value to use. 

441 

442 See Also 

443 -------- 

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

445 """ 

446 # These settings are from the suite2P documentation 

447 TAU_MAP = {'G6s': 1.5, 'G6m': 1., 'G6f': .7, 'default': 1.5} 1 ctx1n

448 _, subject, *_ = session_path_parts(self.session_path) 1 ctx1n

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

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

451 key = match.groups()[0] if match else 'default' 1 ctx1n

452 return TAU_MAP.get(key, TAU_MAP['default']) 1 ctx1n

453 

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

455 """ 

456 Create the ops dictionary for suite2p based on metadata. 

457 

458 Parameters 

459 ---------- 

460 meta: dict 

461 Imaging metadata. 

462 

463 Returns 

464 ------- 

465 dict 

466 Inputs to suite2p run that deviate from default parameters. 

467 """ 

468 

469 # Computing dx and dy 

470 cXY = np.array([fov['Deg']['topLeft'] for fov in meta['FOV']]) 1 ctx1e

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

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

473 

474 # Currently supporting z-stacks but not supporting dual plane / volumetric imaging, assert that this is not the case 

475 if np.any(nXnYnZ[:, 2] > 1): 1 ctx1e

476 raise NotImplementedError('Dual-plane imaging not yet supported, data seems to more than one plane per FOV') 

477 

478 sW = np.sqrt(np.sum((np.array([fov['Deg']['topRight'] for fov in meta['FOV']]) - np.array( 1 ctx1e

479 [fov['Deg']['topLeft'] for fov in meta['FOV']])) ** 2, axis=1)) 

480 sH = np.sqrt(np.sum((np.array([fov['Deg']['bottomLeft'] for fov in meta['FOV']]) - np.array( 1 ctx1e

481 [fov['Deg']['topLeft'] for fov in meta['FOV']])) ** 2, axis=1)) 

482 pixSizeX = nXnYnZ[:, 0] / sW 1 ctx1e

483 pixSizeY = nXnYnZ[:, 1] / sH 1 ctx1e

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

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

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

487 

488 # Computing number of unique z-planes (slices in tiff) 

489 # FIXME this should work if all FOVs are discrete or if all FOVs are continuous, but may not work for combination of both 

490 slice_ids = [fov['slice_id'] for fov in meta['FOV']] 1 ctx1e

491 nplanes = len(set(slice_ids)) 1 ctx1e

492 

493 # Figuring out how many SI Rois we have (one unique ROI may have several FOVs) 

494 # FIXME currently unused 

495 # roiUUIDs = np.array([fov['roiUUID'] for fov in meta['FOV']]) 

496 # nrois = len(np.unique(roiUUIDs)) 

497 

498 db = { 1 ctx1e

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

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

501 'fast_disk': '', # TODO 

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

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

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

505 'keep_movie_raw': False, 

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

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

508 'nimg_init': 400, 

509 'combined': False, 

510 'nonrigid': True, 

511 'maxregshift': 0.05, # default = 1 

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

513 'block_size': [128, 128], 

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

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

516 'mesoscan': True, 

517 'nplanes': nplanes, 

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

519 'nchannels': nchannels, 

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

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

522 'slices': slice_ids, # this tells us which FOV corresponds to which tiff slices 

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

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

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

526 'dx': dx, 

527 'dy': dy 

528 } 

529 

530 return db 1 ctx1e

531 

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

533 """ 

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

535 

536 Parameters 

537 ---------- 

538 run_suite2p: bool 

539 Whether to run suite2p, default is True. 

540 rename_files: bool 

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

542 use_badframes: bool 

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

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

545 

546 Returns 

547 ------- 

548 list of pathlib.Path 

549 All files created by the task. 

550 """ 

551 import suite2p 1 ctx1e

552 

553 """ Metadata and parameters """ 1 ctx1e

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

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

556 collections = sorted(set(f.parts[-2] for f in meta_files)) 1 ctx1e

557 # Check there is exactly 1 meta file per collection 

558 assert len(meta_files) == len(list(self.session_path.glob(self.device_collection))) == len(collections) 1 ctx1e

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

560 if len(rawImagingData) > 1: 1 ctx1e

561 meta = self._check_meta_data(rawImagingData) 

562 else: 

563 meta = rawImagingData[0] 1 ctx1e

564 # Get default ops 

565 ops = suite2p.default_ops() 1 ctx1e

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

567 db = self._create_db(meta) 1 ctx1e

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

569 for k, v in kwargs.items(): 1 ctx1e

570 if k in ops.keys() or k in db.keys(): 1 ctx1e

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

572 db[k] = v 1 ctx1e

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

574 self.kwargs = {**self.kwargs, **db} 1 ctx1e

575 

576 """ Bad frames """ 1 ctx1e

577 # exptQC.mat contains experimenter QC values that may not affect ROI detection (e.g. noises, pauses) 

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

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

580 qc_paths = sorted(map(str, filter(Path.exists, qc_paths))) 1 ctx1e

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

582 if len(exptQC) > 0: 1 ctx1e

583 frameQC, frameQC_names, _ = self._consolidate_exptQC(exptQC) 

584 else: 

585 _logger.warning('No frame QC (exptQC.mat) files found.') 1 ctx1e

586 frameQC = np.array([], dtype='u1') 1 ctx1e

587 frameQC_names = pd.DataFrame(columns=['qc_values', 'qc_labels']) 1 ctx1e

588 

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

590 # badframes.mat contains QC values that do affect ROI detection (e.g. no PMT, lens artefacts) 

591 badframes = np.array([], dtype='i8') 1 ctx1e

592 total_frames = 0 1 ctx1e

593 # Ensure all indices are relative to total cumulative frames 

594 for m, collection in zip(rawImagingData, collections): 1 ctx1e

595 badframes_path = self.session_path.joinpath(collection, 'badframes.mat') 1 ctx1e

596 if badframes_path.exists(): 1 ctx1e

597 raw_mat = loadmat(badframes_path, squeeze_me=True, simplify_cells=True)['badframes'] 

598 badframes = np.r_[badframes, raw_mat + total_frames] 

599 total_frames += m['nFrames'] 1 ctx1e

600 if len(badframes) > 0 and use_badframes is True: 1 ctx1e

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

602 

603 """ Suite2p """ 1 ctx1e

604 # Create alf it is doesn't exist 

605 self.session_path.joinpath('alf').mkdir(exist_ok=True) 1 ctx1e

606 # Remove existing suite2p dir if it exists 

607 suite2p_dir = Path(db['save_path0']).joinpath('suite2p') 1 ctx1e

608 if suite2p_dir.exists(): 1 ctx1e

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

610 # Run suite2p 

611 if run_suite2p: 1 ctx1e

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

613 

614 """ Outputs """ 1 ctx1e

615 # Save and rename other outputs 

616 if rename_files: 1 ctx1e

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

618 else: 

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

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

621 out_files = [f for f in out_files if f.name in [f[0] for f in self.output_files]] 1 ctx1e

622 return out_files 1 ctx1e

623 

624 

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

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

627 

628 priority = 40 (empty)

629 job_size = 'small' (empty)

630 

631 @property (empty)

632 def signature(self): (empty)

633 signature = { 2 ctx1ik

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

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

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

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

638 ('rawImagingData.times_scanImage.npy', self.device_collection, True, True), # register raw 

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

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

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

642 } 

643 return signature 2 ctx1ik

644 

645 def _run(self): (empty)

646 """ 

647 Extract the imaging times for all FOVs. 

648 

649 Returns 

650 ------- 

651 list of pathlib.Path 

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

653 

654 """ 

655 # TODO function to determine nFOVs 

656 try: 2 ctx1ik

657 alf_path = self.session_path / self.sync_collection 2 ctx1ik