diff --git a/src/fox_lab_to_nwb/__init__.py b/src/fox_lab_to_nwb/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/fox_lab_to_nwb/assets/camera_views.png b/src/fox_lab_to_nwb/assets/camera_views.png new file mode 100644 index 0000000..04455a9 Binary files /dev/null and b/src/fox_lab_to_nwb/assets/camera_views.png differ diff --git a/src/fox_lab_to_nwb/assets/trialanalysis_photron.m b/src/fox_lab_to_nwb/assets/trialanalysis_photron.m new file mode 100644 index 0000000..d2ba1aa --- /dev/null +++ b/src/fox_lab_to_nwb/assets/trialanalysis_photron.m @@ -0,0 +1,383 @@ +%INPUT: single trial folder, plot (0 or 1) +%analyzes all data within a file and saves struct as .pcf +%can use this per trial, called as part of "folderanalysis" +%toplot and tosave should =1 if yes + +function trialdata = trialanalysis_photron(trial_path, env, toplot, tosave) + +%if not all the arguments, redefine +%we ain't about checking for special cases +if nargin < 3 + trial_path = uigetdir(); + toplot = 1; + env = 40; %CHANGE ENVELOPE SIZE HERE if running single trial + tosave = 0; +end + +trialdata = struct; + +%load in data +%daq +fd = dir([trial_path filesep '*.fly2']); +ddata = load(fullfile(fd.folder, filesep, fd.name),'-mat'); +ddata = ddata.rec; + +%haltere data +hd = dir([trial_path filesep '*_haltamp.mat']); +hdata = load(fullfile(hd.folder, filesep, hd.name)); + +%antenna/head data +ad = dir([trial_path filesep '*_antennatracking*.csv']); +%make a variable to know if head fixed or head free +% headfree = 0; +% if isempty(ad) +% ad = dir([trial_path filesep '*_headtracking*.csv']); +% headfree = 1; +% end +adata = readtable(fullfile(ad.folder, filesep, ad.name), 'NumHeaderLines',1); + +%wing data +wd = dir([trial_path filesep '*_PROC.mat']); +wdata = load(fullfile(wd.folder, filesep, wd.name)); + +%antenna columns: +%1 = coords (frame number), then x y likelihood for each +%2-4 = L_ant; 5-7 = R_ant; 8-10 = L_base; 11-13 = R_base +%head tracking: +%14-16 = L_OC; 17-19 = R_OC; 20-22 = L_pOR; 23-25 = R_pOR; +%26-28 = L_aVT; 29-31 = R_aVT + +%ASSIGN DATA VARIABLES +%daq +%trial info +duration = ddata.trial.duration; +stimlength = ddata.trial.stimlength; +flynum = ddata.trial.flynum; +trialnum = ddata.trial.flytrial; +volt = ddata.trial.stimvolt; +%data +stamps = ddata.daq.tstamps; +ctrig = ddata.daq.data(:,2); %fastec cam trigger +%otrig = ddata.daq.data(:,3); %opto trigger +ptrig = ddata.daq.data(:,9); %photron trigger +wtrig = ddata.daq.data(:,10); %wind trigger + +%NOT CURRENTLY USING +%including in struct though so it's there if needed? +%otrig = ddata.daq.data(:,3); %opto trigger +dwbf = ddata.daq.data(:,6); +%wbf data from wingbeat amplifier needs to be multiplied by 100: +% 1V = 100Hz +dwbf = dwbf*100; %daq wingbeat freq, added 7/17 +dwbaL = ddata.daq.data(:,4); %temp +% twbaL = twbaL - mean(twbaL(1:2499)); %set to 0 for denoising +dwbaR = ddata.daq.data(:,5); +% twbaR = twbaR - mean(twbaR(1:2499)); +% hutchenL = ddata.daq.data(:,7); +% hutchenR = ddata.daq.data(:,8); + +%ALIGN CAMERAS TO DAQ - CHECK WITH MIKE AGAIN? +%find trigger times +ctrigtime = find(ctrig>3,1)/ddata.daq.fs; +ptrigtime = find(ptrig>3,1)/ddata.daq.fs; + +%topcam (currently skipping sidecam!) +meta_t = fastecMetaReader(fullfile(trial_path, 'TOPCAM_000000.txt')); +ts_t = linspace(1/meta_t.fs, meta_t.numframes/meta_t.fs, meta_t.numframes); +ts_t = ts_t-(ts_t(end)-ctrigtime); +%photron/phantom camera +mp = dir([trial_path filesep 'HALTCAM*.mii']); %since the number changes every time +meta_p = photronMetaReader(fullfile(mp.folder, mp.name)); +ts_p = linspace(1/meta_p.fs, meta_p.numframes/meta_p.fs, meta_p.numframes); +ts_p = ts_p-(ts_p(end)-ptrigtime); + +%% ANTENNAE +%subtract height from y data because y = 0 is at top of image +Lx = adata.L_ant; +Ly = meta_t.height - adata.L_ant_1; +Rx = adata.R_ant; +Ry = meta_t.height - adata.R_ant_1; +%base data just in case +Lxb = adata.L_base; +Lyb = meta_t.height - adata.L_base_1; +Rxb = adata.R_base; +Ryb = meta_t.height - adata.R_base_1; +%find theta +Lt = atan2d((Ly-Lyb),(Lx-Lxb)); +Lt = wrapTo360(Lt); +Rt = atan2d((Ry-Ryb),(Rx-Rxb)); +%find distance between L and R +dant = []; +for j = 1:length(Lx) + dant(j) = pdist2([Lx(j) Ly(j)],[Rx(j) Ry(j)]); +end +dant = dant'; + +%% WING DATA +%denoise daq WBA freq data +dwbf = hampel(dwbf); +%flyalyzer data +wbL = wdata.fly.wingL.angle; +% wbLr = wdata.fly.wingL.root; +wbR = wdata.fly.wingR.angle; +%need to wrap R wing data +wbR = wrapTo180(wbR); +% wbRr = wdata.fly.wingR.root; +%create timestamps: flyalyzer tstamps aren't good, it's always every 20 +wstamps = linspace(ts_t(1),ts_t(end),length(ts_t)/20); + +%zeroing removed 7/17: will be zeroed later! +%zero everything, camera upix is 1351 -> 1350/20 = 67.5 - see +%"cameratimes_phantom.m" for calculation +% wbL = wbL-mean(wbL(1:60)); +% wbR = wbR-mean(wbR(1:60)); + +%7/17 alright we're gonna modify the freq calc bc something's wrong +%NOTE: Mike's code, this sums up image pixels to get the frequency. This +%does include the haltere because it's in the video, but the wings are much +%bigger so it shouldn't make too much difference. +%MODIFIED SUMVID +vr = VideoReader([trial_path filesep 'TOPCAM_000000.avi']); +%create array for video sums +int = zeros(1,vr.NumFrames); +%go through frames +for j = 1:vr.NumFrames + %read image, greyscale but saves as color + img = read(vr,j); + %take only red channel since they're all the same + img = img(:,:,1); + %add up all the pixels in that image and add to time series + int(j)=sum(sum(img)); +end +%pxx is made of FFTs for each timestamp, f is freq, t is time +%zscore gets rid of DC info (low freq) +%window size is 64 w 63 overlap (or 128/127), or use 20/0 so same +%sample number as flyalyzer +%freq range is 150-300, 2000 is fps, 'reassigned' does a thing +%reassigned is kinda like contrast stretching? +[pxx,f,t]=spectrogram(zscore(int),64,63,150:300,2000,'reassigned'); +%time/freq ridge is pretty much the wbf plot, has a .001 penalty to +%make sure there's not much difference between y values (so it doesn't +%jump) +%NOTE: there might be some weird harmonics if the penalty isn't small +%enough, can check by plotting t,calc_wbf +calc_wbf = tfridge(pxx,f,.001); + +%NOTE: there might be some weird harmonics if the penalty isn't small +%enough, can check by plotting t,calc_wbf over top of the spectrogram to +%check if it lines up + % spectrogram(zscore(int),128,127,150:300,2000,'yaxis','reassigned'); + % plot(t,calc_wbf,'r') +%add NaNs back in to beginning so same length as timestamps +calc_wbf = [NaN((length(int) - length(calc_wbf)),1); calc_wbf]; + +%% Haltere Data +%in ampalyzer, haltere base needs to be to the right of the haltere +%base, not in the center! otherwise lots of weird artefacts +%haltere data +%get coordinates and subtract y from height of frame (y=0 at top) +p_height = meta_p.height; +%haltere and wing root +hroot = hdata.haltrpos; +hroot(2) = p_height-hroot(2); +wroot = hdata.wingrpos; +wroot(2) = p_height-wroot(2); +hx = hdata.haltpos(:,1); +hy = p_height - hdata.haltpos(:,2); +%this bit is from Mike and some trial and error... +%find a transformation matrix that turns xy in haltere root to 0,0 and xy in +%wing root to 1,0 +tform = fitgeotrans([hroot;wroot],[0,0;1,0],'nonreflectivesimilarity'); +%use that transformation matrix for the tracked points +[hxn, hyn] = transformPointsForward(tform,hx,hy); +%get angles and wrap +Ht = atan2d(hyn,hxn); %PHOTRON +%get AMPLITUDE (envelope function) +%envelope of 14 looks best for phantom, about 40 for photron -DEFINED AT TOP THIS IS WHY +%WE'RE GOING THROUGH FOLDERS DAMMIT +[Htup,Htlow] = envelope(Ht,env,'peak'); +%find peaks the other way... 10k fps at 200 Hz peak-to-peak should be 50, +%so half a stroke is 25, should be just above that... +[Htu_pks, Htu_locs, Htu_w] = findpeaks(Ht, 'MinPeakHeight',20, 'MinPeakDistance',30); +[Htl_pks, Htl_locs, Htl_w] = findpeaks(-Ht, 'MinPeakHeight',20, 'MinPeakDistance',30); +Htl_pks = -Htl_pks; +% %distance, normalized to distance between wing and root +% %removed 6/28 bc not a great time? kept Hd though for reference +% %alternates between top and bottom of wing stroke; every other peak has a +% %hutchen +Hd = pdist2([0 0],[hxn,hyn]); +% %dips at top of stroke.... +% [Hdup,Hdlow] = envelope(Hd,env,'peak'); +% [Hdu_pks, Hdu_locs, Hdu_w] = findpeaks(Hd, 'MinPeakHeight',.4, 'MinPeakDistance',30); + +%get FREQUENCY +%spectrogram function basically makes an fft (power x function: pxx) for +%each time t, samples every 256 points (255 overlapping) +% 64 & 128 is too small because fps is 10k! +%to be optimal range +%detrend gets rid of all the low freq info (dc info?) +%looking for first set of freq and no harmonics (max wbf is around 230?) +[pxx,f,~]=spectrogram(detrend(Ht),256,255,0:300,meta_p.fs); +%find time frequency ridge (kinda like the avg?) +Hfreq = tfridge(pxx,f,.001); +%since spec uses sliding window, first set of overlap is not used, add NaNs +%to make same size as timestamps +Hfreq = [NaN((length(Ht) - length(Hfreq)),1); Hfreq]; + +if isequal(toplot,1) %plot data + close all + figure('WindowState','maximized'); + t = tiledlayout(2,3); + + nexttile + hold on + %yyaxis left + plot(wstamps,wbR,'r') + %for some reason this is dashed if you don't specify it... + plot(wstamps,wbL,'-b') + % yyaxis right + % plot(stamps,wbf) + title('Wingbeat Amp R (red) & L (blue)') + hold off + + %7/17 look at freq + nexttile + hold on + yyaxis left + ylabel('WB Amp') + plot(stamps,dwbf) + yyaxis right + ylabel('Calc') + plot(ts_t,calc_wbf) + + hold off + nexttile + hold on + plot(ts_t,dant) + title('Distance btwn Antennae') + hold off + + nexttile + hold on + yyaxis left + plot(ts_t,Lt,'b') + ylabel("Ltheta") + yyaxis right + plot(ts_t,Rt,'r') + ylabel("Rtheta") + title('Antennal angle') + hold off + + nexttile + hold on + plot(ts_p,Hfreq) + title('Haltere Freq') + hold off + + % nexttile + % hold on + % yyaxis left + % plot(ts_p,(Htup-Htlow),'r') + % ylim([150 200]) + % yyaxis right + % plot(ts_p,(Hdup-Hdlow),'b') + % ylim([.35 .65]) + % title('Haltere Envelope Change angle (red) & length (blue)') + % hold off + + nexttile + hold on + envelope(Ht,env,'peak') + title('Haltere angle') + hold off + + % nexttile + % hold on + % plot(Ht) + % plot(Htu_locs,Htu_pks,'o') + % plot(Htl_locs,Htl_pks,'o') + % title('Haltere angle peaks') + % hold off + + + % nexttile + % hold on + % envelope(Hd,env,'peak') + % plot(Hdu_locs,Hdu_pks,'o') + % plot(Hdl_locs,Hdl_pks,'o') + % title('Haltere dist') + % hold off + + title(t,{"Fly: " + num2str(flynum) + " Trial: " + num2str(trialnum), ... + "Stim " + num2str(stimlength) + " Voltage: " + num2str(volt), ... + "Envelope: " + num2str(env)}); +end + + %% Add to struct + %note: nothing is normalized to the mean yet! This is only the raw + %data! + %not added: hutchens + %trial info + trialdata.info.duration = duration; + trialdata.info.stimlength = stimlength; + trialdata.info.flynum = flynum; + trialdata.info.trialnum = trialnum; + trialdata.info.envelope = env; + trialdata.info.volt = volt; + %daq data + trialdata.trial.tstamps = stamps; + trialdata.trial.ctrig = ctrig; + %trialdata.trial.otrig = otrig; + trialdata.trial.wbaL = dwbaL; + trialdata.trial.wbaR = dwbaR; + trialdata.trial.wbf = dwbf; + trialdata.trial.ptrig = ptrig; + trialdata.trial.wtrig = wtrig; + %camera data + trialdata.ant.tstamps = ts_t'; + trialdata.ant.dist_ant = dant; + trialdata.ant.Rt = Rt; + trialdata.ant.Lt = Lt; + trialdata.wing.wbL = wbL; + trialdata.wing.wbR = wbR; + trialdata.wing.wstamps = wstamps; + trialdata.wing.wbf = calc_wbf; + trialdata.halt.tstamps = ts_p'; + trialdata.halt.Ht = Ht; + trialdata.halt.Htup = Htup; + trialdata.halt.Htlow = Htlow; + trialdata.halt.Hdist = Hd; + % trialdata.halt.Hdup = Hdup; + % trialdata.halt.Hdlow = Hdlow; + trialdata.halt.Hfreq = Hfreq; + trialdata.halt.Htu_pks = Htu_pks; + trialdata.halt.Htu_locs = Htu_locs; + trialdata.halt.Htu_w = Htu_w; + trialdata.halt.Htl_pks = Htl_pks; + trialdata.halt.Htl_locs = Htl_locs; + trialdata.halt.Htl_w = Htl_w; + % trialdata.halt.Hdu_pks = Hdu_pks; + % trialdata.halt.Hdu_locs = Hdu_locs; + % trialdata.halt.Hdu_w = Hdu_w; + % trialdata.halt.Hdl_pks = Hdl_pks; + % trialdata.halt.Hdl_locs = Hdl_locs; + % trialdata.halt.Hdl_w = Hdl_w; + + % %head data if applicable TO DO + % if headfree == 1 + % trialdata.head.pitch = []; + % trialdata.head.yaw = []; + % trialdata.head.roll = []; + % end + + if tosave == 1 + %% save + %save analysis so it runs faster if you're playing around with the graphs + %date & time + %dt = string(datetime('now'),'yyMMdd_HHmmss'); + %save file name (Tsh file) + flycross = extractBetween(trial_path,'trials\','\'); + flycross = convertCharsToStrings(flycross); + sfn = convertCharsToStrings(trial_path) + filesep + flycross + '_f' + flynum + '_r' + trialnum + '_trial.pcf'; + save(sfn,'trialdata'); + end \ No newline at end of file diff --git a/src/fox_lab_to_nwb/assets/windtrial.m b/src/fox_lab_to_nwb/assets/windtrial.m new file mode 100644 index 0000000..de9bcb5 --- /dev/null +++ b/src/fox_lab_to_nwb/assets/windtrial.m @@ -0,0 +1,181 @@ +function rec = windtrial +%AMS - edited from "windtrial_phantom" in "Antennae_Opto_0623" +%daq setup is for opto trials - may need to modify later +%to use with new photron camera + +%returns struct (rec) as fly2 file + +%need to reset daq, don't need to clear all b/c you're in a function +daqreset + +%% DEFINE VARIABLES HERE FOR EACH TRIAL +%fly info for save file +flynum = 16; +flytrial = 1; +flycross = 'TESTING'; +%directory name UPDATE IF NEEDED +datafolder = "C:\Users\Foxlab\Dropbox\Fox_lab\Amy\Antennae_windtrials"; + +%stimulus info (in seconds) + +stimtime = .5; %length of stimulus +stimint = 1; %interval between stimuli/after stimulus +stimwait = 0.25; %time before first stim +stimrep = 1; %repetitions + +%leaving opto stuff in FOR NOW because of how daq is set up! +%set stimvolt to 0 so it doesn't go off +stimvolt = 0; + +%% END VARIABLE EDITS + + +%% MESSAGE BOX +%this is just copy/pasted from the help page for the inputdlg function +%not as good as a gui like xdrift b/c doesn't save overwritten variables +prompt = {'Cross', 'Fly Number','Trial Number', 'Stimulus Length', 'Inter-Stimulus Interval',... + 'Lag Time', 'Repetitions', 'Opto Volt'}; +dlgtitle = 'Trial Data'; +fieldsize = [1 30; 1 30; 1 30; 1 30; 1 30; 1 30; 1 30; 1 30]; +%default values +definput = {flycross, num2str(flynum), num2str(flytrial), num2str(stimtime), num2str(stimint),... + num2str(stimwait), num2str(stimrep), num2str(stimvolt)}; +%creates a cell array +answer = inputdlg(prompt,dlgtitle,fieldsize,definput); +%if there's something empty then cancel trial +if isempty(answer) + return +end + +%rewrite data +flycross = answer{1}; +flynum = str2double(answer{2}); +flytrial = str2double(answer{3}); +stimtime = str2double(answer{4}); +stimint = str2double(answer{5}); +stimwait = str2double(answer{6}); +stimrep = str2double(answer{7}); +stimvolt = str2double(answer{8}); + + +%define total stim length +stimlength = stimwait + stimrep*(stimtime+stimint); + +%% Make DAQ objects & channels +%INITIATE DAQ +d = daq('ni'); +d.Rate = 10000; %10k rate + +%ADD DAQ CHANNELS +fastectrig = addoutput(d, "Dev2","ao0","Voltage"); +optotrig = addoutput(d,"Dev2","ao1","Voltage"); %DISCONNECTED +phantomtrig = addoutput(d,"Dev2","ao2","Voltage"); +windtrig = addoutput(d, "Dev2","ao3","Voltage"); %WIND + +CamSync = addinput(d, "Dev2", "ai0","Voltage"); +CamTrigger = addinput(d, "Dev2", "ai1","Voltage"); +OptoTrigger = addinput(d, "Dev2", "ai2","Voltage"); +LWingbeatAmp = addinput(d, "Dev2", "ai3","Voltage"); +RWingbeatAmp = addinput(d, "Dev2", "ai4","Voltage"); +WingbeatFreq = addinput(d, "Dev2", "ai5","Voltage"); +LHutchen = addinput(d, "Dev2", "ai6","Voltage"); +RHutchen = addinput(d, "Dev2", "ai7","Voltage"); +%there's no AI 8-15 on the daq? +PTrigger = addinput(d, "Dev2", "ai16","Voltage"); +WTrigger = addinput(d, "Dev2", "ai17","Voltage"); + +%START SUB-FUNCTION (this would return "scandata" and not rec!) + +%CREATE TRIGGER/INPUT VECTORS +%camera trigger is the stimlength multiplied by the rate, then UP about +%.25sec after? could use the stimtime for the trigger but if you want the +%stimulus to be longer then the file will be too big. +%all input vectors need to be single columns, and all the same size! +%camera trigger - goes up a the end of the trial +%multiply by 3.3 because the fastec cameras want a voltage of 3.3 +trigtime = 0.25; %time after trigger/when trig is 1! using 250 ms so +camtrig = [zeros(stimlength*d.Rate,1); ones(trigtime*d.Rate,1)]*3.3; +%Phantom camera is TTL! needs 5V to trigger, won't trigger at 3.3 +phantomtrig = [zeros(stimlength*d.Rate,1); ones(trigtime*d.Rate,1)]*5; +totallength = stimlength + trigtime; +%KEEPING THIS IN FOR NOW, OPTO LIGHT NOT IN PLACE +%opto trig, might have variable sizes and lengths so add with a for loop +otrig = [zeros(stimwait*d.Rate,1)]; %add rate time +for i=1:stimrep %COULD EDIT + otrig = [otrig; ones(stimtime*d.Rate,1); zeros(stimint*d.Rate,1)]; +end +%need to add to the end so it's the same length as the trigger +otrig = [otrig; zeros(0.25*d.Rate,1)]; +%multiply by voltage if we want to vary intensity (change LED driver to +%"mod") +otrig = otrig*stimvolt; + +%WIND TRIGGER +wtrig = [zeros(stimwait*d.Rate,1)]; %add rate time +for i=1:stimrep %COULD EDIT + wtrig = [wtrig; ones(stimtime*d.Rate,1); zeros(stimint*d.Rate,1)]; +end +wtrig = [wtrig; zeros(0.25*d.Rate,1)]; +windvolt = 9; %solenoid meant to run at 7V, doesn't react until 9V +%solenoid power source should be 12V btw +wtrig = wtrig*windvolt; + + +%indicator LED + +%% GO! +%create output matrix, needs to be MxN matrix, where M is num data scans +%and N is num output channels +outdata = [camtrig otrig phantomtrig wtrig]; +%readwrite function: inScanData = readwrite(d,outScanData); +scandata = readwrite(d,outdata,"OutputFormat","Matrix"); + +%% CREATE STRUCT AND SAVE FILE +%create struct, needs to follow xdrift naming! +rec = struct; +%daq data +rec.daq.fs = d.Rate; +rec.daq.data = scandata; +rec.daq.tstamps = linspace(0,totallength,totallength*d.Rate)'; +rec.daq.channelnames = ["CamSync" "CamTrigger" "OptoTrigger" "LWingBeatAmp" "RWingBeatAmp" "WingBeatFreq" "LHutchen" "RHutchen" "PTrigger" "WindTrigger"]; + +%need to create base file name here for save +%need to convert flycross from char vector to string +flycross = convertCharsToStrings(flycross); +ffn = flycross + '_' + string(datetime('now'),'yyMMdd_HHmmss') + '_f' + string(flynum) + '_r' + string(flytrial); + +%trial info, can also be figured out from daq data +rec.trial.info = ffn; %trial info (filename) +rec.trial.duration = stimlength; +rec.trial.stimlength = stimtime; +rec.trial.interval = stimint; +rec.trial.delay = stimwait; +rec.trial.reps = stimrep; +rec.trial.flynum = flynum; +rec.trial.flytrial = flytrial; +rec.trial.stimvolt = stimvolt; +rec.trial.windvolt = windvolt; + +%SAVE +%create trial folder, this will also create the cross folder if it doesn't +%exist! just make sure the datafolder is correct! +%directory file name +dfn = datafolder + filesep + flycross + filesep + ffn; +mkdir(dfn); +%save file name (fly2 file) +sfn = dfn + filesep + ffn + '.fly2'; +save(sfn,'rec'); + +%END SUB-FUNCTION HERE? + +%PLOT TO CHECK WBA +wbaL = rec.daq.data(:,4); +wbaR = rec.daq.data(:,5); +wbf = rec.daq.data(:,6); +plot(wbaL - wbaL(1)) +hold on +plot(wbaR - wbaR(1)) +plot(wbf - wbf(1)) +hold off + +end diff --git a/src/fox_lab_to_nwb/behavior.py b/src/fox_lab_to_nwb/behavior.py new file mode 100644 index 0000000..bf973de --- /dev/null +++ b/src/fox_lab_to_nwb/behavior.py @@ -0,0 +1,107 @@ +from pathlib import Path +from typing import Optional + +from neuroconv import BaseDataInterface +from pydantic import FilePath +from pynwb import NWBFile +from pynwb import TimeSeries +from pymatreader import read_mat + + +class BehaviorInterface(BaseDataInterface): + def __init__(self, file_path: FilePath): + self.file_path = Path(file_path) + + def add_to_nwbfile( + self, + nwbfile: NWBFile, + metadata: Optional[dict] = None, + ) -> None: + + mat = read_mat(self.file_path) + recording_structure = mat["rec"] + daq_struct = recording_structure["daq"] + + timestamps = daq_struct["tstamps"] + + # Extracted from the matlab manually + # neither scipy or pymatreader can read string characters from the .mat file + channel_names = [ + "CamSync", + "CamTrigger", + "OptoTrigger", + "LWingBeatAmp", + "RWingBeatAmp", + "WingBeatFreq", + "LHutchen", + "RHutchen", + "PTrigger", + ] + + # TODO: figure how how to store this + # Synchronization signals + cam_sync = daq_struct["data"][:, 0] + cam_trigger = daq_struct["data"][:, 1] + opto_trigger = daq_struct["data"][:, 2] + ptrigger = daq_struct["data"][:, 8] + + # Behavior signals + left_wing_beat_amplitude = daq_struct["data"][:, 3] + right_wing_beat_amplitude = daq_struct["data"][:, 4] + wing_beat_frequency = daq_struct["data"][:, 5] + + unit = "tbd" + description = "tbd" + left_wing_beat_amplitude_time_series = TimeSeries( + name="LeftWingBeatAmplitudeTimeSeries", + data=left_wing_beat_amplitude, + unit=unit, + timestamps=timestamps, + description=description, + ) + + description = "tbd" + right_wing_beat_amplitude_time_series = TimeSeries( + name="RightWingBeatAmplitudeTimeSeries", + data=right_wing_beat_amplitude, + unit=unit, + timestamps=timestamps, + description=description, + ) + + description = "tbd" + unit = "Hz" # TODO: Figure this out, the values in the plot are around 3 but should be higher for flies + wing_beat_frequency_time_series = TimeSeries( + name="WingBeatFrequencyTimeSeries", + data=wing_beat_frequency, + unit=unit, + timestamps=timestamps, + description=description, + ) + + nwbfile.add_acquisition(left_wing_beat_amplitude_time_series) + nwbfile.add_acquisition(right_wing_beat_amplitude_time_series) + nwbfile.add_acquisition(wing_beat_frequency_time_series) + + # TODO: Ask Ben if this nesting makes sense? probably not + # time_series = [left_wing_beat_amplitude_time_series, right_wing_beat_amplitude_time_series, wing_beat_frequency_time_series] + # behavioral_time_series_container = BehavioralTimeSeries(name="BehavioralTimeSeries", time_series=time_series) + # nwbfile.add_acquisition(behavioral_time_series_container) + + # Not clear what are those signals, haltere? + lhutchen = daq_struct["data"][:, 6] + rhutchen = daq_struct["data"][:, 7] + + unit = "tbd" + description = "tbd" + lhutchen_time_series = TimeSeries( + name="LHutchenTimeSeries", data=lhutchen, unit=unit, timestamps=timestamps, description=description + ) + + description = "tbd" + rhutchen_time_series = TimeSeries( + name="RHutchenTimeSeries", data=rhutchen, unit=unit, timestamps=timestamps, description=description + ) + + nwbfile.add_acquisition(lhutchen_time_series) + nwbfile.add_acquisition(rhutchen_time_series) diff --git a/src/fox_lab_to_nwb/conversion_notes.md b/src/fox_lab_to_nwb/conversion_notes.md new file mode 100644 index 0000000..d857987 --- /dev/null +++ b/src/fox_lab_to_nwb/conversion_notes.md @@ -0,0 +1,172 @@ + +# General comments +Point-person: Amy Streets (axs2909@case.edu) and Kris Lea (kxl786@case.edu) + +https://sites.google.com/site/cwrufoxlab/people + +Paper is still not out, but according to the authors it should be similar to this one: + +https://www.biorxiv.org/content/10.1101/2024.03.13.583703v1 + +The data is available here: + + + +## Trial structure + +Folder name: +- Tshx18D07_240124_115923_f3_r1 + +Info here: +Tshx18D07 = cross +240124 = date +115923 = time +f3 = fly number +r1 = trial repeat number + + +Organized with folder hierarchy, with one folder for all trials of a single type/cross containing individual trial folders (ex. \Tshx18D07\Tshx18D07_240124_115923_f3_r1 where the + +# DeepLabCut output files +config files missing but Video and DLC output files are available + +### Pickle (metadata) + +Here is an example output: + +``` +{'data': {'start': 1710281222.2714453, + 'stop': 1710281243.7364442, + 'run_duration': 21.464998960494995, + 'Scorer': 'DLC_resnet50_antennatrackingMar11shuffle1_100000', + 'DLC-model-config file': {'stride': 8.0, + 'weigh_part_predictions': False, + 'weigh_negatives': False, + 'fg_fraction': 0.25, + 'mean_pixel': [123.68, 116.779, 103.939], + 'shuffle': True, + 'snapshot_prefix': 'C:\\Users\\FoxLab-Whiskey\\Desktop\\antennatracking-AMS-2024-03-11\\dlc-models\\iteration-0\\antennatrackingMar11-trainset95shuffle1\\test\\snapshot', + 'log_dir': 'log', + 'global_scale': 0.8, + 'location_refinement': True, + 'locref_stdev': 7.2801, + 'locref_loss_weight': 1.0, + 'locref_huber_loss': True, + 'optimizer': 'sgd', + 'intermediate_supervision': False, + 'intermediate_supervision_layer': 12, + 'regularize': False, + 'weight_decay': 0.0001, + 'crop_pad': 0, + 'scoremap_dir': 'test', + 'batch_size': 8, + 'dataset_type': 'imgaug', + 'deterministic': False, + 'mirror': False, + 'pairwise_huber_loss': True, + 'weigh_only_present_joints': False, + 'partaffinityfield_predict': False, + 'pairwise_predict': False, + 'all_joints': [[0], [1], [2], [3]], + 'all_joints_names': ['L_ant', 'R_ant', 'L_base', 'R_base'], + 'dataset': 'training-datasets\\iteration-0\\UnaugmentedDataSet_antennatrackingMar11\\antennatracking_AMS95shuffle1.mat', + 'init_weights': 'C:/Users/FoxLab-Whiskey/Desktop/antennatracking-AMS-2024-03-11\\dlc-models\\iteration-0\\antennatrackingMar11-trainset95shuffle1\\train\\snapshot-100000', + 'net_type': 'resnet_50', + 'num_joints': 4, + 'num_outputs': 1}, + 'fps': 30.0, + 'batch_size': 8, + 'frame_dimensions': (480, 640), + 'nframes': 4351, + 'iteration (active-learning)': 0, + 'training set fraction': 0.95, + 'cropping': False, + 'cropping_parameters': [0, 640, 0, 480]}} +``` + +# Static TIFF files for confocal imaging +Not yet available + +# Synchronization signal with Spike2 +Not yet available + +# Intracellular electrophysiology data including the voltage trace and stimulus trace + +Not yet available + +# Files on data shared + +### Tshx18D07_240124_115923_f3_r1.fly2 +This is supposed to be a struct with DAQ data and the wingbeat analysis. Can't open it with matlab online unless extension is changed to .mat + + Channel names: +- CamSync +- CamTrigger +- OptoTrigger +- LWingBeatAmp +- RWingBeatAmp +- WingBeatFreq +- LHutchen +- RHutchen +- PTrigger +- WindTrigger + +These can't be extracted with pyton from matlab because they are strings. Are they always the same order? + +## Cameras + +### Fastec + +Files: +* TOPCAM_000000.avi +* SIDE_000000.avi + +Format is avi, metadata is in xml + +https://www.fastecimaging.com/ + +TopCam has a DLC analysis for the following body parts: + +bodyparts: L_ant, R_ant, L_base, R_base, + +Sidecam does not have DLC analysis + +### Phantom +Files: +* XZ_1_186.mp4 + +Format is mp4, metadata is in xml + +This camera has an associated DLC analysis for the following body parts: + +bodyparts: haltere + +### Photron + +Files: +* Tshx18D07_240124_115923_f3_r1_down20.mp4 [This seems to be the same as topcam] + + +Not sure about this. The `.mii` metadata is not available and there is this extra file: + + +It also does not have a corresponding DLC analysis + +### Graphical example for the trials +![example](./assets/camera_views.png) + +See the example, we have SideCam, TopCam and XZ_1_186.mp4 + +The other vide is Tshx18D07_240124_115923_f3_r1_down20.mp4 which I guess was downsampled. + + +## `Tshx18D07_f3_r1_trial.Tsh` +This should be trial data. +Unable to read in matlab online unless extension is changed to .mat + +## Tshx18D07_240124_115923_f3_r1_down20_PROC.mat +This has the analysis for the wingbeat redone. + + +## Tshx18D07_ant_top_f3_r1_500ms_dvProject.mat +Not described \ No newline at end of file diff --git a/src/fox_lab_to_nwb/metadata.yaml b/src/fox_lab_to_nwb/metadata.yaml new file mode 100644 index 0000000..1e20f17 --- /dev/null +++ b/src/fox_lab_to_nwb/metadata.yaml @@ -0,0 +1,21 @@ +NWBFile: + keywords: + - Fly + - Wingbeat + - Behavior + related_publications: + - https://doi.org/great_publication or link to APA or MLA citation of the publication + session_description: + A rich text description of the experiment. Can also just be the abstract of the publication. + experiment_description: 'Task: Rapid serial visual presentation (RSVP).' + institution: Case Western Reserve University + lab: Fox + experimenter: + - Streets, Amy + - Lea, Kriss + - Fox, Jessica + surgery: TBD +Subject: + sex: F + species: Drosophila melanogaster + description: fly diff --git a/src/fox_lab_to_nwb/streets_conversion.py b/src/fox_lab_to_nwb/streets_conversion.py new file mode 100644 index 0000000..e552ba2 --- /dev/null +++ b/src/fox_lab_to_nwb/streets_conversion.py @@ -0,0 +1,130 @@ +from pathlib import Path +import time +from datetime import datetime +from typing import Optional +from zoneinfo import ZoneInfo + +import parse +from neuroconv.utils import dict_deep_update, load_dict_from_file +from neuroconv import ConverterPipe +from neuroconv.datainterfaces import DeepLabCutInterface, VideoInterface + +from fox_lab_to_nwb.behavior import BehaviorInterface + + +def run_trial_conversion(trial_folder_path: Path, output_dir_path: Optional[Path] = None, verbose: bool = True): + + if verbose: + start_time = time.time() + + if output_dir_path is None: + output_dir_path = Path.home() / "conversion_nwb" + + # Parse metadta from the trial_folder_path name + # Define the correct parsing pattern + pattern = r"{cross}_{date}_{time}_f{fly_number}_r{trial_repeat_number}" + parsed_metadata = parse.parse(pattern, trial_folder_path.name) + + trial_repeat_number = parsed_metadata["trial_repeat_number"] + cross = parsed_metadata["cross"] # TODO, where in subject metadata should this be? Example Tshx18D07 + fly_number = parsed_metadata["fly_number"] + session_date = parsed_metadata["date"] + session_time = parsed_metadata["time"] + + subject = f"Fly{fly_number}" + + session_id = f"{subject}_{session_date}_{session_time}" + nwbfile_path = output_dir_path / f"{session_id}.nwb" + + # Behavior interface + format = "fly2" # Authors said in email this might change + daq_file_name = f"{cross}_{session_date}_{session_time}_f{fly_number}_r{trial_repeat_number}.{format}" + file_path = trial_folder_path / daq_file_name + behavior_interface = BehaviorInterface(file_path=file_path) + + # Video Interface + # TODO: are the names of the video files always the same per trial? we need + # More trials to find out + + side_cam_name = "SideCam_000000.avi" + file_path = trial_folder_path / side_cam_name + side_cam_interface = VideoInterface(file_paths=[file_path], metadata_key_name="SideCam") + + top_camera_name = "TOPCAM_000000.avi" + file_path = trial_folder_path / top_camera_name + + top_cam_interface = VideoInterface(file_paths=[file_path], metadata_key_name="TopCam") + + haltere_camera_name = "XZ_1_186.mp4" + + file_path = trial_folder_path / haltere_camera_name + haltere_cam_interface = VideoInterface(file_paths=[file_path], metadata_key_name="BackCam") + + # DLC interface + top_cam_dlc_file_name = "TOPCAM_000000DLC_resnet50_antennatrackingMar11shuffle1_100000.h5" + top_cam_file_path = trial_folder_path / top_cam_dlc_file_name + top_cam_dlc_interface = DeepLabCutInterface(file_path=top_cam_file_path) + + haltere_cam_dlc_file_name = "XZ_1_186DLC_resnet50_haltereMar13shuffle1_100000.h5" + haltere_cam_file_path = trial_folder_path / haltere_cam_dlc_file_name + haltere_cam_dlc_interface = DeepLabCutInterface(file_path=haltere_cam_file_path) + + data_interface = { + "Behavior": behavior_interface, + "SideCam": side_cam_interface, + "TopCam": top_cam_interface, + "BackCam": haltere_cam_interface, + "DeepLabCutTopCam": top_cam_dlc_interface, + "DeepLabCutHaltereCam": haltere_cam_dlc_interface, + } + + converter = ConverterPipe(data_interfaces=data_interface) + + # Add datetime to conversion + metadata = converter.get_metadata() + session_datetime = datetime.strptime(f"{session_date}_{session_time}", "%y%m%d_%H%M%S") + session_start_time = session_datetime.astimezone(ZoneInfo("America/New_York")) + metadata["NWBFile"]["session_start_time"] = session_start_time + + # Update default metadata with the editable in the corresponding yaml file + editable_metadata_path = Path(__file__).parent / "metadata.yaml" + editable_metadata = load_dict_from_file(editable_metadata_path) + metadata = dict_deep_update(metadata, editable_metadata) + + subject_metadata = metadata["Subject"] + subject = "subject" + subject_metadata["subject_id"] = f"{subject}" + + # Run conversion, this adds the basic data to the NWBFile + conversion_options = { + "DeepLabCutTopCam": {"container_name": "PoseEstimationTopCam"}, + "DeepLabCutHaltereCam": {"container_name": "PoseEstimationHaltereCam"}, + } + + converter.run_conversion( + metadata=metadata, + nwbfile_path=nwbfile_path, + conversion_options=conversion_options, + overwrite=True, + ) + + if verbose: + stop_time = time.time() + conversion_time_seconds = stop_time - start_time + if conversion_time_seconds <= 60 * 3: + print(f"Conversion took {conversion_time_seconds:.2f} seconds") + elif conversion_time_seconds <= 60 * 60: + print(f"Conversion took {conversion_time_seconds / 60:.2f} minutes") + else: + print(f"Conversion took {conversion_time_seconds / 60 / 60:.2f} hours") + + +if __name__ == "__main__": + + verbose = True + data_path = Path("/home/heberto/cohen_project/Sample data/Fox Lab") + assert data_path.exists(), f"Folder {data_path} does not exist" + trial_folder_path = data_path / "Tshx18D07_240124_115923_f3_r1" + assert trial_folder_path.exists(), f"Folder {trial_folder_path} does not exist" + + run_trial_conversion(trial_folder_path=trial_folder_path, verbose=verbose)