%% Processing after Suite2p analysis
% Written by Seung-Eon Roh (soh3@jhmi.edu), last updated September 26, 2020.
% This code is based on the data generated by Suite2p Python version.

% From GUI one can save the whole signal data set by selecting File-'Save
% to mat file (*.mat)'.
clear all; close all;
load('Fall.mat'); F=F-Fneu % Subtract background signals
csvwrite('F_minus_neu.csv',F); % Save it


%% Basic signal processing and plotting

% Selected ROIs
A=[0 3]; % enter ROI numbers that you are interested and saw in the GUI.
A=A+1; F_select=F(A,:)'; S_select=spks(A,:)';
csvwrite('F_select.csv', F_select); csvwrite('S_select.csv', S_select);

% Z-score Normalization
% F signals and S(spike data) signals are normalized.
[r c] = size(F_select); F_std=zscore(F_select); S_std=zscore(S_select);
F_sort=sort(F_std); F_base=mean(F_sort(1:round(r*0.1),:));% Lower 10% average of F
var_corr=F_base-min(S_std(:,:)); S_corr=(S_std-abs(var_corr))*0.5;
S_corr=S_corr+abs(min(S_corr)); S_corr(S_corr < 0.5) = 0;% Exclude values under 0.5
csvwrite('F_std_sel.csv', F_std); csvwrite('S_std_sel.csv', S_std); % Save normalized data

% Plot settings
[r, c, z] = size(F_select); Hz=32; % put your own sampling rate for 'Hz'
sec_length = r/Hz;
xval = 1/Hz:1/Hz:sec_length; m = 5; % 'm' determines how far from other signal data

% F_std display with separation of 'm'
f_display = zeros(r,c,z);
    for i = 1:z
        for j = 1:c
            f_display(:,j,i) = F_std(:,j,i) -m*(j-1);
        end
    end
    for i = 1:z
        figure
        plot(xval, f_display(:,:,i))
        title(i); xlabel('Time (sec)');
    end
    
% Generate raster plot from S_corr
vec=(1:1:r); S_corr_raster=zeros(r,c);
for i = 1 : c
S_corr1=S_corr(:,i);
[pks,locs] = findpeaks(S_corr1);
vec_exc=setdiff(vec,locs');
S_corr1(vec_exc)=0;
S_corr_raster(:,i)=S_corr1;
end
csvwrite('S_raster.csv', S_corr_raster);

% S_corr_raster display (upon the F_std plot above)
f_display = zeros(r,c,z);
    for i = 1:z
        for j = 1:c
            f_display(:,j,i) = S_corr_raster(:,j,i) -m*(j-1);
        end
    end
    for i = 1:z
        hold
        %f_display(f_display==0)=nan;
        plot(xval, f_display(:,:,i),'color',[0, 0, 0]);
        title(i); xlabel('Time (sec)');
    end
    
    
       
%% Spontaneous response amplitudes and traces
% (1) This is to obtain amplitudes for all spontaneous responses and collect
% signal data around the event. Apuff responses excluded.

% Compute time-range for APuff to be excluded for analysis
% Below it computes air puff every 5 sec from at 10 sec of recording.
AP_no=(r-Hz*10)/Hz/5; locs_AP=[(5*1+5)*Hz-2:(5*1+5)*Hz+3];
for i=2:AP_no;
 locs_AP= [locs_AP (5*i+5)*Hz-2:(5*i+5)*Hz+3];
end
locs_AP=sort(locs_AP); %

% Extract spontaneous amplitudes
amp_spon=zeros(r,c);
for i =1:c;
  casignal=S_corr(:,i);
  [pks locs]=findpeaks(casignal);
  locs=setdiff(locs,locs_AP); % excludes timepoints of APuff
  amp_spon(locs,i)=S_corr(locs,i);
end

% Take 1st events within 20 frames to avoid resurge of signals
% This is needed because multiple adjacent events compound the
% interpretation.
amp_spon_1stpk=amp_spon;
for j = 1:c;
    for i = 1:r-9;
      if amp_spon(i,j)>0.1;
       amp_spon_1stpk(i+1:i+19,j)=0;
       end;
    end; 
end
amp_spon_1stpk(1:5,:)=0; % delete events during first 5 frames
amp_spon_1stpk(r-30:r,:)=0; % delete events during last 30 frames; prevents time-point mismatch

% Amplitude option 1
amp_no_spon=sum(amp_spon_1stpk~=0);amp_mean=sum(amp_spon_1stpk) ./ amp_no_spon;
amp_spon_total=amp_spon_1stpk; amp_spon_total(amp_spon_total==0) = [];
% 'amp_spon_total' contains all amplitudes of spontaneous events.

% (2) extract spon F_signal 
% This is to collect all signal traces from before the number of (pre) to
% after the number of (post)
pre=4; post=22; total=pre+post+1; % set length of the responses to be collected
sig_spon_no=nnz(amp_spon_1stpk); sig_spon=zeros(total,sig_spon_no); % generate empty matrix

[row col]=find(amp_spon_1stpk(:,:)); sig_spon=zeros(total,numel(row));
ind_spon=[row col]; % compute index for spon responses
for i=1:numel(row);
    sig_spon(1:total,i)=F_std(ind_spon(i,1)-pre:ind_spon(i,1)+post,ind_spon(i,2));
end % Collect all traces of spontaneous responses

sig_spon_mean=mean(sig_spon,2); figure; plot(sig_spon); % plot(sig_spon_mean);
save('sig_spon', 'sig_spon'); csvwrite('sig_spon.csv', sig_spon); % save traces

% Amplitude option 2
% Below is amplitudes collected directly from the traces. Please choose
% this option if you want to compare spon vs Apuff responses. Otherwise
% Amplitude option 1 is enough.
amp_spon_total=max(sig_spon(1:10,:)); amp_spon_mean=mean(amp_spon_total);
save('amp_spon_total', 'amp_spon_total'); csvwrite('sig_spon.csv', sig_spon);
csvwrite('amp_spon.csv', amp_spon_total);



%% Apuff response amplitudes and traces based on APuff templete

% IMPORTANT- Generate tempelet for APuff
temp_AP=ones(c,AP_no); % Please execute this to generate templete and then MANUALLY switch 1 to 0, where the Apuff responses failed as you see the traces.
save('tempelet_AP','temp_AP'); % You can save the templete (optional).
% load('tempelet_AP');

% Definition of Apuff times
locs_AP_evt=[];
for i=1:AP_no;
    locs_AP_evt=[locs_AP_evt; (5*i+5)*Hz:(5*i+5)*Hz+5];
end

% exporting traces for air puff times
APtraces=[]; twdw=5;
for j=1:AP_no
    APt=locs_AP_evt(j,1);
    tmp=[];
    for k=1:c
        if temp_AP(k,j)==1
            [m,id]=max(F_std(APt:APt+twdw-1,k)); % Find time position of maximum value
            MAXt=APt+id-1;
            tmp=[tmp F_std(MAXt-5:MAXt+twdw+16,k)]; % Shift trace to set maximum position = initial position
        else
            tmp=[tmp NaN(twdw+22,1)]; % Unless validated in temp_AP, fill the column with NaN
        end
    end
    APtraces=[APtraces tmp];
end
figure; plot(APtraces); % plot the traces for Apuff
csvwrite('sig_AP.csv', APtraces); save('sig_AP', 'APtraces'); % save it
amp_AP=max(APtraces(1:10,:)); amp_AP = amp_AP(~isnan(amp_AP));
% amp_AP contains all amplitudes of Apuff responses
save('amp_AP', 'amp_AP'); csvwrite('amp_AP.csv', amp_AP);
amp_AP_mean=nanmean(amp_AP); 



%% Plot spon vs AP traces
sig_avg_spon=mean(sig_spon'); sig_avg_AP=nanmean(APtraces');
figure; plot(sig_avg_spon); hold on; plot(sig_avg_AP);
amp_spon_AP=[amp_spon_mean amp_AP_mean]; figure; bar(amp_spon_AP);
