workflow.m

%% declare folders with data

files{1}.VSD_dir = '/storage/laur/Data_2/CoattiAlex/Experiments_2015/L17_210115/RedShirtImaging';
files{4,1}.VSD_dir = '/storage/laur/Data_2/CoattiAlex/Experiments_2015/turtle_6March/Ale6March_SampleD';
files{8,1}.VSD_dir = '/storage/laur/Data_2/CoattiAlex/Experiments_2015/turtle_6March/Ale6March_SampleG';
% For PC
files{1}.VSD_dir_PC = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\L17_210115\RedShirtImaging';
%% declare ROIs and/or frames to read out from data

% index = {[78:87],[70:75],[1:30000];[78:87],[149:154],[1:30000];...
%     [128:137],[110:115],[1:30000];[177:186],[72:77],[1:30000];...
%     [177:186],[151:157],[1:30000];}
% index = {[72:159],[84:146],[1:16000]};

index = {[1:256],[1:256],[1:30000]};

querymode = 'fullframe';  %'rois'; %'fullframe';
meanframe = 1; % saves mean F per frame for full recording and timings of led on for synch
write_to_disk = 1; % if 1 then deinterleave based on vertically flipping one half of the camera fov
darkframe = 1; % if 1 uses last frame to substract a dark frame from all frames
exp_id = 1;
rep_id = 2;
PC_flag = 1; % 1 when working on microsoft PC

files = read_Tsm_redshirt (files,exp_id,rep_id,index,querymode,meanframe,write_to_disk,darkframe,PC_flag);
% version 2: work around problem of reading tsm header
files = read_Tsm_redshirt2 (files,exp_id,rep_id,index,querymode,meanframe,write_to_disk,darkframe,PC_flag);


import matlab.io.*;
fptr = fits.openFile('GabazineVSD2001.tsm');

%% Read in data from ANDOR recordings

recording_dir = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\Q18_260315\Andor_VSD';
cd(recording_dir); content_folder = dir;

% specify folder name where data of interested are saved
get_name = 'sample_4';
folder_ids = find(cell2mat(cellfun(@(x) x==1, {content_folder.isdir},'UniformOutput',0))); 
select_folder = find(~cellfun('isempty', strfind({content_folder(folder_ids).name},get_name)));
cd(content_folder(folder_ids(select_folder)).name); list = ls;
for l = 1:length(list)
    file_flag(l) = ~isempty(strfind(list(l,:),'sif'));
end
filenames = list(find(file_flag),:);
% for all input parameters that sifread_ale can accept see help sifread_ale
% optional key parameter: 
% 1) frames_to_read_startime = vector of frame numbers from which the function starts reading data
% 2) window_read = number of frames read in from corresponding frame_startime
% notes: the function is set to default read_mode (1) therefore takes in
% the optional inputs {frames_to_read_startime,window_read}; if these are
% not declared than by default all frames will be read in using read_mode=1

for i = [1 8 9 12 13 14] %1:length(filenames)
    absfilepath = strcat(recording_dir,filesep,get_name,filesep,filenames(i,:));
    [newdata, meta] = sifread_ale(absfilepath);
    data_rec(i).data = newdata;
    data_rec(i).metadata = meta;
    clear newdata meta
end

%% Average fluorescence over time
exp_id = [1 8 9 12 13 14]; nframes = zeros(1,length(exp_id)); c = 0;
for i = exp_id
    c = c+1;
    nframes(c) = data_rec(i).metadata.NumberImages;
    exp_name{c} = data_rec(i).metadata.FileName;
end
avg_F_dt = NaN(length(exp_id),max(nframes));
for i = 1:length(exp_id) 
    tempdata = data_rec(exp_id(i)).data;
    avg_F_dt(i,1:size(tempdata,3)) = mean(mean(tempdata,1),2);
end
figure, plot(avg_F_dt'), legend(exp_name)
%% Convert ANDOR sif file into tiff



%% Determine frames for synching with MEA & onset/offset illumination

exp_ids = 1;
rep_ids = 1:5;
shift = 100; % an abrupt change is defined when a differen in intensity 
    % between two successive frames exceeds the the mean fluctuation observed 
    % between successive frames
centremass_option = 0; % if 1 then once events where an abrupt change in illumination 
    % is found it gives a time of the event based on the centre of mass of the
    % intensities of light

files = find_transients(files,exp_ids,rep_ids,shift,centremass_option)

%% Align MEA and VSD recordings

% Create object for MCRack Data
folder_data = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\Q18_260315\MC_Rack'
get_name = 'sample_4';
cd(folder_data)
% specify folder name where data of interested are saved
folder_ids = find(cell2mat(cellfun(@(x) x==1, {content_folder.isdir},'UniformOutput',0))); 
select_folder = find(~cellfun('isempty', strfind({content_folder(folder_ids).name},get_name)));
cd(content_folder(folder_ids(select_folder)).name); list = ls; file_flag = [];
for l = 1:length(list)
    file_flag(l) = ~isempty(strfind(list(l,:),'mcd'));
end
filenames = list(find(file_flag),:);

i= 8;
MC = MCRackRecording(strtrim(filenames(i,:)));

[Dig_data,trig_ms] = getDigitalData(MC);
Digital_data = squeeze(Dig_data(1,1,:));
Digital_timestamps = trig_ms;

vsd_expid = 12;
% camera parameters
camera_acquisition_frequency = 1/data_rec(vsd_expid).metadata.KineticCycleTime;
movie_duration_ms = data_rec(vsd_expid).metadata.NumberImages/camera_acquisition_frequency * 1000;
n_frames = data_rec(vsd_expid).metadata.NumberImages;
% mea parameters
mea_acquisition_frequency = MC.samplingFrequency;
first_high = find(Digital_data,1,'first');
last_high = first_high + find( Digital_data(first_high:first_high+movie_duration_ms* mea_acquisition_frequency/1000) , 1,'last');
all_high = numel(find(Digital_data(first_high:last_high)));
step = all_high/n_frames;

MEA_trace_high = squeeze(squeeze(LFP_channel(:,:,find(Digital_data(first_high:last_high)))));
MEA_trace_scaled = ( MEA_trace_high - mean(MEA_trace_high) ) ./ max(abs(MEA_trace_high - mean(MEA_trace_high)));
vsd_timepoints = [1:step:all_high];
vsd_trace_scaled = ( avg_F_dt(4,:) - mean(avg_F_dt(4,:)) ) ./ max(abs( avg_F_dt(4,:) - mean(avg_F_dt(4,:)) ));

% find high-low transitions in digital data (very inefficient)
p=find(dig_data==1);
p=p(2.035e6:end);


[LFP_channel t_ms] = MC.getData(236,Digital_timestamps(first_high),Digital_timestamps(last_high)-Digital_timestamps(first_high));

[LFP_channel t_ms] = MC.getData(1,3.622e5,(4.162-3.622)*1e5);

%% Mask for VSD recording
i=1;
avg_image = mean(data_rec(i).data,3);
threshold_intensity = 2000;
th_mask = avg_image;
th_mask(find(avg_image < threshold_intensity)) = 0;

[r,c] = find(logical(th_mask));

width = data_rec(1).metadata.width;
height = data_rec(1).metadata.height;
pxl_ind = reshape([1:width*height],width,height);

reduced_fov = pxl_ind(min(r)-1:max(r),min(c):max(c));
bw_cropped_mask = logical(th_mask_cropped);

% to be handled over to VSD_rise_delay
p_r = 5; p_c = 3;
step_r = 1:p_r:size(reduced_fov,1); step_c = 1:p_c:size(reduced_fov,2);
%% visualize delays in VSD signals 

exp_id = 1; rep_id = 1:5;
p_r = 8;
p_c = 8;
darkframe = 1; % if 1 uses trace{1,2}, if 0 uses trace{1,1} and substracts mean(frames(1:100))

files = VSD_rise_delay(files,exp_id,rep_id,p_r,p_c,index,darkframe);

for j = 1:size(delays,2) 
    delaytemp = delays(:,j);
    delay_plot(:,:,j) = reshape(delaytemp,length(step_r),length(step_c));
end

for j = 1:size(delays,2)
figure, imagesc(delay_plot(:,:,j)), colorbar
axis square
end

   %% compare traces from different trials (create a function)
   
interval{1} = [707:6706,9711:15706];
interval{2} = [527:527+5997,9526:9526+5997];
interval{3} = [490:490+5997,9489:9489+5997];
interval{4} = [657:657+5997,9657:9657+5997];
interval{5} = [674:674+5997,9672:9672+5997];
interval{6} = [489:489+5997,9489:9489+5997];
interval{7} = [524:524+5997,9523:9523+5997];


rois = {[70:80],[40:50];[30:40],[40:50];[78:88],[10:20];[55:65],[25:35]}

for i = 8
    traces = NaN(length(interval{1}),size(rois,1),7);
    F0 = NaN(length(interval{1}),size(rois,1),7);
    dFoF = NaN(length(interval{1}),size(rois,1),7);
    for j = 1:7
    
        temp = files{i,j}.RawData.traces{1,1};
        DF = mean(temp(:,:,1:300),3);
        movie = temp - repmat(int16(DF),1,1,size(temp,3));
%         figure, plot(squeeze(mean(mean(movie,1),2)))
        movie_chopped = movie(:,:,interval{j});
        W = 1500;
        for r = 1:size(rois,1)
            traces(:,r,j) = mean(mean(movie_chopped(rois{r,1},rois{r,2},:),1),2);
            tracetemp = squeeze(mean(mean(movie_chopped(rois{r,1},rois{r,2},:),1),2));
            tot_length = size(tracetemp,1);
            xic = tracetemp(1:floor(W/2));
            xfc = tracetemp(tot_length - floor(W/2)+1:tot_length);
            F0temp = fastmedfilt1d(tracetemp, W, xic, xfc);
            F0(:,r,j) = F0temp;
            dFoF(:,r,j) = (tracetemp - F0temp) ./ F0temp;
        end
    end
end