function [tSRD tSRDfull] = AnalyzeGCaMP_2017(tSRDfull)
%Function takes in an n x 5 matrix where the first column is time (s), second
%column is size of object (in square microns), 3rd column is mean ratio channel,
%4th column is X position of centroid (in microns). The 5th column is Y position
% of centroid (in microns). Sixth column in output file will be delta R over R once
%it is calculated.
%
%If no argument is provided, will prompt user to choose a comma delineated (.csv)
%file or folder of .csv files. Each file should have 5 columns and no words or
%headers.
%tSRD stands for time size ratio and delta R over R. tSRDfull(:,3) will
%be adjusted by a rolling average of size aveSize (default rolling
%average = 1; e.g. no rolling average. The resulting matrix is tSRD.
%If there is no input, grab a file, otherwise set tSRDfull to input
if nargin == 0
[name,path,~] = uigetfile('.csv','multiselect','on');
if ~isequal(class(name),'cell')
name = {name};
end
tSRDfull = cell(length(name),1);
for i = 1:length(name)
tSRDfull{i} = csvread(strcat(path,char(name{i})));
end
end
%Create variable to hold tSRDfull and to manipulate
temptSRD = tSRDfull;
%Create the size of tSRD
tSRD = cell(length(temptSRD),1);
%Size of rolling average
aveSize = 1;
%baseline in decimal. This will be used to calculate the lowest X
%percent of data points. If a recording has a lot of Ca2+ activity,
%this value should be set lower to avoid negative ?R/R values.
cutOff = 0.1;
%makes figure to show progress of peak finding and ?R/R
figure;
%Do the following for every file. loop through each file
for z = 1:length(temptSRD)
%weight temptSRD by multiplying area column with ratio and XY
%centroid columns. This will form an 'area-adjusted ratio' and
%area-adjusted centroid position that will allow larger objects to
%count for more of the relative ratio and centroi position compared
%to smaller objects. This area adjustment will be removed at a later step
temptSRD{z}(:,3) = temptSRD{z}(:,2).*temptSRD{z}(:,3);
temptSRD{z}(:,4) = temptSRD{z}(:,2).*temptSRD{z}(:,4);
temptSRD{z}(:,5) = temptSRD{z}(:,2).*temptSRD{z}(:,5);
%loop through and eliminate doubled time points by adding together the size
%of the two objects and averaging the signal between them. Loop
%will continue until doubled timepoints are consolidated
for i = length( temptSRD{z}(:,1) ) : -1 : 2
if temptSRD{z}(i,1) == temptSRD{z}(i-1,1)
temptSRD{z}(i-1,2) = temptSRD{z}(i-1,2) + temptSRD{z}(i,2);
temptSRD{z}(i-1,3) = temptSRD{z}(i-1,3) + temptSRD{z}(i,3);
temptSRD{z}(i-1,4) = temptSRD{z}(i-1,4) + temptSRD{z}(i,4);
temptSRD{z}(i-1,5) = temptSRD{z}(i-1,5) + temptSRD{z}(i,5);
temptSRD{z}(i,:) = [];
end
end
%go through the ratio, X centroid, and Y centroid values, dividing
%by the cumulative area of the now-summed objects. Reverses
%area-adjustment incorporated earlier duing object consolidation
temptSRD{z}(:,3) = temptSRD{z}(:,3)./temptSRD{z}(:,2);
temptSRD{z}(:,4) = temptSRD{z}(:,4)./temptSRD{z}(:,2);
temptSRD{z}(:,5) = temptSRD{z}(:,5)./temptSRD{z}(:,2);
%Set up tSRD as a subset of temptSRD and add an extra column of zeros
tSRD{z} = [temptSRD{z}(1:end-aveSize+1,:) zeros(length(temptSRD{z})-aveSize+1,1)];
%turn tSRD's area column to a rolling average
for k = 1:length(tSRD{z}(:,1))
tSRD{z}(k,2) = mean(temptSRD{z}(k:k+aveSize-1,2));
end
%turn tSRD's ratio column to a rolling average
for j = 1:length(tSRD{z}(:,1))
tSRD{z}(j,3) = mean(temptSRD{z}(j:j+aveSize-1,3));
end
%turn tSRD's X column to a rolling average
for k = 1:length(tSRD{z}(:,1))
tSRD{z}(k,4) = mean(temptSRD{z}(k:k+aveSize-1,4));
end
%turn tSRD's Y column to a rolling average
for k = 1:length(tSRD{z}(:,1))
tSRD{z}(k,5) = mean(temptSRD{z}(k:k+aveSize-1,5));
end
%sort tSRD to find baseline
sorttSRD = sort(tSRD{z}(:,3));
%find the index at which the cutoff would sit
cutOffIndex = floor(cutOff*length(tSRD{z}));
%calculate baseline as mean of the bottom cutOff of values
baseline = mean(sorttSRD(1:cutOffIndex));
%calculate delta R over R
tSRD{z}(:,6) = (tSRD{z}(:,3) - baseline)/baseline;
%write the matrix to a file
if ~exist(strcat(path,'/analyzed'),'dir')
mkdir(path,'analyzed')
end
%Writes 6-column output to a new folder called 'analyzed', appending
%'a-' to original filename. Data output columns are: Time, Area, Ratio,
%X centroid, Y centroid, and finally ?R/R of recording.
%Output has additional precision for very long recordings where the
%hundreths place of the time column would get rounded. This is not
%necessary if timepoint (e.g. frame of the recording) is used instead.
dlmwrite(strcat(path,'/analyzed/a-',char(name(z))),tSRD{z},'precision','%6.3f')
%Using findpeaks algorithm to find transients in column 6 (?R/R).
%Returns the amplitude (pks), timepoint (locs), width (w), and
%prominance (like amplitude but based on local area, not raw
%intensity). pks (?R/R); timepoint (s), width (s), prominance
%(?R/R). MinPeakDistance is the minimum amount of time,
%in seconds, between peaks. MinPeakDistance = 1 means peaks will only
%be detected if they are 1s or more apart. It elminates peaks with two or more
%maxima with the same value, although it doesn't fix the width or
%other calculated features of those peaks. MinPeakProminence is in
%units of delta-R/R and measure the fraction of the local baseline
%rather than raw ratio from baseline. Writes a matrix with the prefix 'peaks'
%using PeakFile with that information.
%Note, the top iteration of the findpeaks algorithm writes the
%peak matrix, the bottom iteration draws the graphs with annotated
%peaks. You must change the paramters of both for them to agree.
[pks,locs,w,p] = findpeaks(tSRD{z}(:,6),tSRD{z}(:,1),'Annotate','extents','WidthReference','halfheight','MinPeakDistance',1,'MinPeakProminence',0.20);
PeakFile{z} = [locs,pks,w,p];
%write the peak matrix to a file including times, amplitudes,
%widths and prominances.
if ~exist(strcat(path,'/peaks'),'dir')
mkdir(path,'peaks')
end
dlmwrite(strcat(path,'/peaks/peaks-',char(name(z))),PeakFile{z},'precision','%6.3f')
%write the deltaR/R traces to a file
if ~exist(strcat(path,'/traces'),'dir')
mkdir(path,'traces')
end
%plot and save un-annotated deltaR/R in eps format readable by
%Illustrator, etc.
plot(tSRD{z}(:,1),tSRD{z}(:,6));
saveas(gcf,strcat(path,'/traces/plain-',char(name(z)),'.eps'),'epsc');
%plot and save annotated deltaR/R in eps format readable by
%Illustrator, etc.
findpeaks(tSRD{z}(:,6),tSRD{z}(:,1),'Annotate','extents','WidthReference','halfheight','MinPeakDistance',1,'MinPeakProminence',0.20);
saveas(gcf,strcat(path,'/traces/anno-',char(name(z)),'.eps'),'epsc');
end
end