normandCalciumspikeanalyse.m

%-- 28-01-10 10:35 --%

%Program interpolating a fast calcium response
%Data traces should be organised in couloumns. Sampling time should be in the first coloumn%
%First row should be the name of the ROIs or stimulations%
%data should be in a .txt format and placed in a directory folder of the
%same name.
%The starttime.xls document should contain a spread sheet with the names of
%the experiment folder in the first column, the start of stimulation time
%in the second column and the stimulation frequency in the third column.

clear all

%Parameter that changes from experiment to experiment%

stimlength =15;
Directory = 'Q:\TransNeuro\barbara\fra lasaf\e221210\';
Expdate = 'e221210';
ROInumber= 20; %number of ROIs, should match the number of coloumns in your imported data. 
resamplefrek=0.01; %The time intervals between interpolated data points.
Method = 'linear'; %Interpolation method
Scanspeed = 0.085;


%%
%getting stimulationlist
[inflist, stimlist, inforaw] =xlsread([Directory Expdate 'starttime.xlsx']);
Values = zeros(size(stimlist,1),ROInumber+1,8);
Roiname = cell(1,ROInumber+1);
ENtext={'ROIs'};
HzC= cell(1,1);


for e=1:size(stimlist,1);%For each stimulation
ExpName =cell2mat(stimlist(e,1));
tstart =inflist(e,1);%baseline period from beginning of scan to start of repeated stimulation 
StimHz =inflist(e,2);%stimulation frequency

%getting data
A=importdata([Directory ExpName '\' ExpName '.txt']);
data=A.data;
textdata=A.textdata;

%%


%Finding Ca2+ signal spikes%

StimFreq=1/StimHz; %number of seconds between stimulations%
m = round(StimFreq/Scanspeed); %number of scans between stimulations%
numbofroi= (size(data,2))-1;

%%

%Normalising data%

tiden = data (:,1); %coloumn with sampling times%
baseline = data (tiden<tstart,:); %Baseline part of data%
average = mean (baseline); %calculating mean fluorescence in baseline period%
stdev = std (baseline);  %Standard deviationen in fluorescence in baseline period%


i = data (:,2:end); %all data except sampling time%

x = zeros (size(data));%Matrix for normalized data%
stdx = zeros (size(data));%Matrix for normalized data%

for i = 2:size (data,2); 
x(:,i)= ((data (:,i)- average (1,i))/ average (1,i))*100; %Normalising to deltaF/F in percent%
stdx(:,i)= (data (:,i)- average (1,i))/ stdev (1,i); %Matrix to standard deviations of baseline activity%
end

x (:,1) = data (:,1); % sampling times added%
stdx (:,1) = data (:,1); % sampling times added%

% Computed features converted to cell format%
xcell = num2cell(x);
stdxcell = num2cell(stdx);

y = [textdata; xcell]; %adding ROI names to the normalised data%

xlswrite ([Directory 'Stim' ExpName '.xlsx'], y, 'PROdata'); % Saving normalised data%

%%
%Plotting and saving each timeseries as tif and adobe illustrator%
figh = figure(1);
for i = 2:size (x,2); 
plot (x(:,1),x(:,i)); 
xlabel('Time [s]');
ylabel('Ca++ change [%]');
ROInr = textdata{2,i};
saveas (figh, [[Directory ExpName] '\PROstim' ROInr '.tif']);
saveas (figh, [[Directory ExpName] '\PROstim' ROInr '.ai']);
end

%%

%signal spikes calculations%

numbstim = ceil(stimlength/StimFreq); %number of stimulation in total%
tpointer = cell(numbstim, size(x,2)); %prelocating matrix%
spikeresp = cell (numbstim, (size(x,2))+1, 2); %prelocating matrix%
rela_spikeresp = cell (size(spikeresp));
samle = zeros (((StimFreq/0.01)+1), size(x,2));
stdsamle = zeros (((StimFreq/0.01)+1), size(x,2));

%data after each stimulation pooles%

for i = 2:size (x,2); 
    for t=1:numbstim; %For each stimulation at a time%
        
        a = tstart - StimFreq; %defines were in the timeseries signals must be estimated from%
        
        tpointer{t,i}=find(x(:,1)<(a+(t*StimFreq)+StimFreq) & x(:,1)>=(a+(t*StimFreq))); %limits of every spike signal%
         
        spikeresp{t,i,1}=x(tpointer{t,i},i); %Data from ROI coloumns from tpointer
        spikeresp{t,i,2}=stdx(tpointer{t,i},i); %Data from ROI coloumns from tpointer
        spikeresp{t,1,1}=x(tpointer{t,i},1); 
        spikeresp{t,2+numbofroi,1}=x(tpointer{t,i},1)-(StimFreq*t)-a; 
             
        % Removing slow trends in data
        
        if t==1
            addf = tpointer{1,i}(1,1);
            baselinepkt = x (addf-1,i);
            stdbaselinepkt = stdx (addf-1,i);
        else
            addr = tpointer{t-1,i}(end);
            baselinepkt = x(addr,i);
            stdbaselinepkt = stdx(addr,i);
        end      

        rela_spikeresp {t,i,1} = spikeresp{t,i,1}- baselinepkt;
        rela_spikeresp {t,i,2} = spikeresp{t,i,2}- stdbaselinepkt;
        
        rela_spikeresp{t,1,1}=x(tpointer{t,i},1);   %addding sampling time values%
        rela_spikeresp{t,2+numbofroi,1}=x(tpointer{t,i},1)-(StimFreq*t)-a; 
        rela_spikeresp{t,1,2}=x(tpointer{t,i},1); %addding sampling time values%
        rela_spikeresp{t,2+numbofroi,2}=x(tpointer{t,i},1)-(StimFreq*t)-a; %Normalised sampling time values%
        
    end
    
end

relaspike_mat = cell2mat (rela_spikeresp(:,:,1)); 
stdrelaspike_mat = cell2mat (rela_spikeresp(:,:,2)); 
samle2 = zeros ((size(relaspike_mat,1)),(size(x,2)));
stdsamle2 = zeros ((size(relaspike_mat,1)),(size(x,2)));
tid2 = zeros ((size(relaspike_mat,1)),(size(x,2)));
stdtid2 = zeros ((size(relaspike_mat,1)),(size(x,2)));

%Interpolating calcium spikes%

for i = 2:size (x,2); 
    
    ROInr = textdata{2,i};
    ts=timeseries(relaspike_mat(:,i),relaspike_mat(:,2+numbofroi),'Name',ROInr);%timeseries for each ROI with all data%
    ts = setinterpmethod (ts,Method); %defining interpolation method%
    res_ts= resample(ts,0:resamplefrek:StimFreq); %Interpolation%    
    samle (:,i) = res_ts.data;
    samle2 (:,i) = ts.data;
    tid2 (:,i) = ts.time; %getting timevalues from ts%
    stdts=timeseries(stdrelaspike_mat(:,i),stdrelaspike_mat(:,2+numbofroi),'Name',ROInr);%timeseries for each ROI with all data%
    stdts = setinterpmethod (stdts,Method);  %defining interpolation method%
    stdres_ts= resample(stdts,0:resamplefrek:StimFreq); %Interpolation%   
    stdsamle (:,i) = stdres_ts.data;
    stdsamle2 (:,i) = stdts.data;
    stdtid2 (:,i) = stdts.time; %getting timevalues from ts%
    
end

%Setting new timevalues%

samle2 (:,1) = tid2 (:,2);
stdsamle2 (:,1) = tid2 (:,2);

for xx = 1:size (samle,1);
    
    samle (xx,1) = -1*resamplefrek + resamplefrek*xx;
    stdsamle (xx,1) = -1*resamplefrek + resamplefrek*xx;
    
end

samle_cell = num2cell(samle);
samlealt = [textdata; samle_cell]; %adding ROI names%

samle2_cell = num2cell(samle2);
samle2alt = [textdata; samle2_cell];


xlswrite ([Directory 'Spike' ExpName '.xlsx'], samlealt, 'PROspikedata'); % saving the interpolated data%
xlswrite ([Directory 'Spike' ExpName '.xlsx'], samle2alt, 'PROspikerawdata');% saving the interpolated data%

disp ([ExpName 'saving results in' Directory 'PROSpike' ExpName '.xlsx'])

save ([Directory 'PROSpike' ExpName '.mat'],'spikeresp','rela_spikeresp','samle')

disp ([ExpName 'saving results in' Directory 'PROSpike' ExpName '.mat'])

%%
%Plotting and saving each spike as tif and adobe illustrator%

figh = figure(2);

for i = 2:size (samle,2);  
    plot (samle(:,1),samle(:,i),'LineWidth',2); 
    hold on;
    plot (samle2(:,1),samle2(:,i),'r.','MarkerSize',8);
    hold off;
    xlabel('Time [s]');
    ylabel('Ca++ change [%]');
    ROInr = textdata{2,i};
    saveas (figh, [[Directory ExpName] '\PROspike' ROInr '.tif']);
    saveas (figh, [[Directory ExpName] '\PROspike' ROInr '.ai']);
end

disp (['End' ExpName])

end
disp ('The End')