Add files via upload

+helper/bonf_holm.m

@@ -0,0 +1,113 @@
+
+% Bonferroni-Holm (1979) correction for multiple comparisons.  This is a
+% sequentially rejective version of the simple Bonferroni correction for multiple
+% comparisons and strongly controls the family-wise error rate at level alpha.
+%
+% It works as follows:
+% 1) All p-values are sorted in order of smallest to largest. m is the
+%    number p-values.
+% 2) If the 1st p-value is greater than or equal to alpha/m, the procedure
+%    is stopped and no p-values are significant.  Otherwise, go on.
+% 3) The 1st p-value is declared significant and now the second p-value is
+%    compared to alpha/(m-1). If the 2nd p-value is greater than or equal 
+%    to alpha/(m-1), the procedure is stopped and no further p-values are 
+%    significant.  Otherwise, go on. 
+% 4) Et cetera.
+%
+% As stated by Holm (1979) "Except in trivial non-interesting cases the 
+% sequentially rejective Bonferroni test has strictly larger probability of
+% rejecting false hypotheses and thus it ought to replace the classical 
+% Bonferroni test at all instants where the latter usually is applied."
+%
+%
+% function [corrected_p, h]=bonf_holm(pvalues,alpha)
+%
+% Required Inputs:
+%  pvalues - A vector or matrix of p-values. If pvalues is a matrix, it can
+%            be of any dimensionality (e.g., 2D, 3D, etc...).
+%
+% Optional Input:
+%  alpha   - The desired family-wise alpha level (i.e., the probability of
+%            rejecting one of more null hypotheses when all null hypotheses are
+%            really true). {default: 0.05}
+%
+% Output:
+%  corrected_p  - Bonferroni-Holm adjusted p-values. Any adjusted p-values
+%                 less than alpha are significant (i.e., that null hypothesis 
+%                 is rejected).  The adjusted value of the smallest p-value
+%                 is p*m.  The ith smallest adjusted p-value is the max of
+%                 p(i)*(m-i+1) or adj_p(i-1).  Note, corrected p-values can
+%                 be greater than 1.
+%  h            - A binary vector or matrix of the same dimensionality as
+%                 pvalues.  If the ith element of h is 1, then the ith p-value
+%                 of pvalues is significant.  If the ith element of h is 0, then
+%                 the ith p-value of pvalues is NOT significant.
+%
+% Example:
+% >>p=[.56 .22 .001 .04 .01]; %five hypothetical p-values
+% >>[cor_p, h]=bonf_holm(p,.05)
+% cor_p =
+%    0.5600    0.4400    0.0050    0.1200    0.0400
+% h =
+%     0     0     1     0     1
+% 
+% Conclusion: the third and fifth p-values are significant, but not the
+% remaining three.
+%
+% Reference:
+% Holm, S. (1979) A simple sequentially rejective multiple test procedure.
+% Scandinavian Journal of Statistics. 6, 65-70.
+%
+%
+% For a review on contemporary techniques for correcting for multiple
+% comparisons that are often more powerful than Bonferroni-Holm see:
+%
+%   Groppe, D.M., Urbach, T.P., & Kutas, M. (2011) Mass univariate analysis 
+% of event-related brain potentials/fields I: A critical tutorial review. 
+% Psychophysiology, 48(12) pp. 1711-1725, DOI: 10.1111/j.1469-8986.2011.01273.x 
+% http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/mass_uni_preprint1.pdf
+%
+%
+% Author:
+% David M. Groppe
+% Kutaslab
+% Dept. of Cognitive Science
+% University of California, San Diego
+% March 24, 2010
+function [corrected_p, h]=bonf_holm(pvalues,alpha)
+if nargin<1
+    error('You need to provide a vector or matrix of p-values.');
+else
+    if ~isempty(find(pvalues<0,1)),
+        error('Some p-values are less than 0.');
+    elseif ~isempty(find(pvalues>1,1)),
+        fprintf('WARNING: Some uncorrected p-values are greater than 1.\n');
+    end
+end
+if nargin<2
+    alpha=.05;
+elseif alpha<=0
+    error('Alpha must be greater than 0.');
+elseif alpha>=1
+    error('Alpha must be less than 1.');
+end
+s=size(pvalues);
+if isvector(pvalues)
+    if size(pvalues,1)>1
+       pvalues=pvalues'; 
+    end
+    [sorted_p sort_ids]=sort(pvalues);    
+else
+    [sorted_p sort_ids]=sort(reshape(pvalues,1,prod(s)));
+end
+[dummy, unsort_ids]=sort(sort_ids); %indices to return sorted_p to pvalues order
+sorted_p(isnan(sorted_p)) = [];
+m=length(sorted_p); %number of tests
+mult_fac=m:-1:1;
+cor_p_sorted=sorted_p.*mult_fac;
+cor_p_sorted(2:m)=max([cor_p_sorted(1:m-1); cor_p_sorted(2:m)]); %Bonferroni-Holm adjusted p-value
+cor_p_sorted = [cor_p_sorted NaN(1,length(sort_ids)-length(cor_p_sorted))];
+corrected_p=cor_p_sorted(unsort_ids);
+corrected_p=reshape(corrected_p,s);
+h=corrected_p<alpha;
+end

+helper/calculateSNR.m

@@ -0,0 +1,63 @@
+
+
+function SNR  = calculateSNR(sData)
+
+     % peakSNR=Δpeak/ 2*σn
+     % Δpeak is the difference between the biggest spike value 
+     % and σn is the noise SD calculated from nonactive intervals of traces.
+    DFF = sData.imdata.roiSignals(2).dff;
+    [noROIs,noSample] = size(DFF);
+
+    % Low-pass filter fluorescence signals to find noise and actual signal
+
+    DFF_lowpass = NaN(noROIs,noSample);
+
+%     % check frequency distribution of single cells: 
+%     roiNO = 10;
+%     s  = DFF(roiNO,:);
+%     L = length(s);                                      % Vector Length
+%     Ts = 1/31;                                          % Sampling Interval (sec)
+%     Fs = 1/Ts;                                          % Sampling Frequency (Hz)
+%     Fn = Fs/2;                                          % Nyquist Frequency (Hz)
+%     t = 1:Ts:L*Ts;                                      % Time Vector
+%     FTs = fft(s-mean(s))/L;                             % Fourier Transform (Subtract d-c Offset)
+%     Fv = linspace(0, 1, fix(L/2)+1)*Fn;                 % Frequency Vector
+%     Iv = 1:length(Fv);                                  % Index Vector
+% 
+%     figure()
+%     plot(Fv, abs(FTs(:,Iv))*2)
+%     hold on
+%     hold off
+%     grid
+%     axis([0 1 ylim])
+%     
+%     
+    % filter paramters:
+    fpass = 0.1;
+    fs    = 31;
+    for i = 1:noROIs
+        DFF_lowpass(i,:) = lowpass(double(DFF(i,:)),fpass,fs);% filter(Hd,DFF(i,:)); 
+    end
+    noise = DFF - DFF_lowpass;  % subtract the lowpassed data from baseline-corrected signal to get the noise
+
+
+
+%   check noise, DFF and DFF_lowpasse:
+%     figure()
+%     plot(noise(roiNO,:))
+%     hold on
+%     plot(DFF(roiNO,:))
+%     plot(DFF_lowpass(roiNO,:),'k')
+%    
+
+    % calculate SNR 
+    noiseStd = NaN(noROIs,1); % whole session std for each ROI
+    ROIPeak  = NaN(noROIs,1); % the highest peaks in that ROI 
+    SNR      = NaN(noROIs,1); % signal to noise ratio for each ROI
+    for i = 1:noROIs
+        noiseStd(i) = 2*std(noise(i,:),[],2); % can be adjusted to 1*std(PreNoise(i,:),[],2);
+        ROIPeak(i)  = prctile(DFF(i,:),99); % can be adjusted to 95%
+        SNR(i)      = ROIPeak(i) / noiseStd(i);
+    end
+end
+

+helper/count_cells.m

@@ -0,0 +1,13 @@
+load_all_data()
+
+
+for r =     1:length(area)
+    count_area(r) = 0;
+    for l = 1:length(file{r})
+
+       sData = load(fullfile(load_sData_dirct,area{r},'/',file{r}{l},'/final_sData.mat'));
+       count_area(r) = count_area(r)+size(sData.imdata.roiSignals(2).dff,1);
+
+    end
+
+end

+helper/createRMaps.m

@@ -0,0 +1,25 @@
+rmaps = struct();
+
+for k =1:length(data)
+
+    rmaps(k,j).dff  = [];
+    rmaps(k,j).deconv  = [];
+
+    for j = 1:length(data(k).sessionIDs)
+        sData  = mData(k,j).sData;        
+        signal_dff = sData.imdata.roiSignals(2).dff;
+        signal_deconv = sData.imdata.roiSignals(2).deconv;
+
+        for f = 1:size(signal_dff,1)                
+                signalMatr                  = helper.splitInTrialsPos(signal_dff, sData, f);                
+                rmaps(k,j).dff(f,:)        =  nanmean(signalMatr,1);
+
+                signalMatr             = helper.splitInTrialsPos(signal_deconv, sData, f);
+                [signalMatr, ~]        = smoothdata(signalMatr, 2,'gaussian',8);
+
+                rmaps(k,j).deconv(f,:) = nanmean(signalMatr,1);
+
+        end
+
+    end
+end

+helper/extractLickFrequency.m

@@ -0,0 +1,39 @@
+
+
+function [lickFrequencyDs ,lickFrequency]= extractLickFrequency(sData,samplingRate,lickLengthThr,lickFreqThr)
+
+lickSignal          = sData.daqdata.lickSignal;
+lickStartInd        = find(diff(lickSignal) == 1)+1;
+lickEndInd          = find(diff(lickSignal) == -1);
+
+% Correct for licks not fully captured at the beginning or end of the recording
+if lickSignal(1) == 1
+    lickEndInd = lickEndInd(2:numel(lickEndInd));
+end
+if lickSignal(numel(lickSignal)) == 1
+    lickStartInd= lickStartInd(1:numel(lickStartInd)-1);
+end
+
+% convert length in sample points to length in ms
+lickLength  = (lickEndInd - lickStartInd)/(samplingRate/1000); 
+
+% convert to lick frequency
+diffLickStart   = diff(lickStartInd);
+lickFreq        = samplingRate./diffLickStart;
+
+% remove too short lick events (and licks with too high frequency)
+shortLickInd = lickStartInd(find(lickLength< lickLengthThr));
+tooFastLickInd = lickStartInd(find(lickFreq > lickFreqThr));
+
+lickErrorsInd = union(shortLickInd,tooFastLickInd);
+lickStartInd  = setdiff(lickStartInd,lickErrorsInd);
+
+lickEvents(1:numel(lickSignal)) = zeros;
+lickEvents(lickStartInd) = 1;
+
+lickFrequency        = helper.ifreq(lickEvents,samplingRate);
+lickFrequencyDs      = lickFrequency(sData.daqdata.frameIndex(1):323:end);
+
+
+% lickFrequencyDs =[];
+end

+helper/getStats.m

@@ -0,0 +1,58 @@
+function roiStat = getStats(sData, CH)
+%getRoiActivityStats Return struct of different roi activity statistics
+%
+%   roiStat = getRoiActivityStats(sData) Returns struct with statistics
+%   about the roi signal for all rois in sData. 
+%
+%   roiStat = getRoiActivityStats(sData, CH) calculates the roiStat on the
+%   specified channel number CH. CH is an integer. The default value is 2.
+%
+%   Output, roiStat contains the following fields:
+%       peakDff         : Peak dff of the roi time series. 
+%       signalToNoise   : Signal to noise 
+%       activityLevel   : Fraction of time where activity is above noise
+%   
+
+% %     pstr = getFilePath(sData.sessionID, 'roiStat');
+% %     if exist(pstr, 'file')
+% %         roiStat = loaddata(sData.sessionID, 'roiStat');
+% %         return
+% %     end
+
+    if nargin < 2
+        CH = 2; % Most of us are imaging on ch2?
+    end
+
+    dff = double(squeeze(sData.imdata.roiSignals(CH).dff));
+    [nRois, nSamples] = size(dff);
+
+
+    % Get max DFF of all rois.
+    peakDff = max(dff, [], 2);
+
+
+    % Get SNR of all Rois.
+    signalToNoise = zeros(nRois, 1);
+    noiseLevel = zeros(nRois, 1);
+    for i = 1:nRois
+        noiseLevel(i) = real(GetSn(dff(i, :)));
+        signalToNoise(i) = snr(dff(i, :), ones(1,nSamples) * noiseLevel(i));
+    end
+
+   signalToNoise = helper.calculateSNR(sData);
+
+    % Get fraction of time above noise level
+    %dffSmooth = okada(dff, 2); % Smooth with okada filter
+    dffSmooth = smoothdata(dff, 2, 'movmean', 7); % Smooth again with movmean
+
+    isHigh = dffSmooth > noiseLevel;
+
+    activityLevel = sum(isHigh, 2) ./ nSamples;
+
+    roiStat = struct('peakDff', peakDff, ...
+                     'signalToNoise', signalToNoise, ....
+                     'activityLevel', activityLevel);
+
+%     savedata(sData.sessionID, struct('roiStat', roiStat))
+
+end

+helper/ifreq.m

@@ -0,0 +1,16 @@
+function iFrequency = ifreq(signal,fs)
+
+% returns the instantaneous frequency of a discrete signal (containing 0 and 1) with fs sampling rate
+
+iFrequency = nan(size(signal));
+if size(signal,1) > size(signal,2)
+    signal = signal';
+end
+events = find(signal == 1);
+freques = fs./diff(events);
+freques = [freques 0];
+iFrequency(signal == 1) = freques;
+iFrequency(1) = 0;
+iFrequency = fillmissing(iFrequency,'previous'); 
+
+end

+helper/splitInTrialsPos.m

@@ -0,0 +1,13 @@
+function binnedTrialedSingleROI = splitInTrialsPos(signal, sData,idx)
+
+for i = 1:numel(sData.trials.trialStart)-1
+    isInTrial   = sData.trials.trialLength == i;
+    trialSignal = signal(idx,isInTrial);
+
+    for jj = 1:78
+        idxInTrialBin                = sData.behavior.wheelPosDsBinned(isInTrial) == jj;
+        binnedTrialedSingleROI(i,jj) = nanmean(trialSignal(idxInTrialBin));
+    end
+end
+
+end

+plot/createHeatMap.m

@@ -0,0 +1,23 @@
+
+function matrix = createHeatMap(lap, pos, behavior)
+
+    % make all the same length
+    minLength = min([length(lap),length(pos), length(behavior)]);
+    lap = lap(1:minLength);
+    pos = pos(1:minLength);
+    behavior = behavior(1:minLength);
+
+    heatmap = zeros(max(lap)+2,ceil(max(pos))); 
+    occupancymap = zeros(max(lap)+2,ceil(max(pos)));
+    offset = min(lap);
+    for x = 1:length(lap)%pos)
+        x_pos = round(pos(x))+1;
+        if x_pos > size(heatmap,2); x_pos = size(heatmap,2); end
+        y_pos = lap(x)+1+(-offset);
+        heatmap(y_pos,x_pos) = heatmap(y_pos,x_pos) + behavior(x);
+        occupancymap(y_pos,x_pos) = occupancymap(y_pos,x_pos) + 1;
+    end
+
+    matrix = heatmap./occupancymap;
+
+end