iPSCquantv2.m

function [] = iPSCquantv2
    %This function loads a list dialogue in your current directory and
    %loads the given sweeps as a structure using the wavesurfer library
    %ws.loadDataFile. Wavesurfer (https://github.com/JaneliaSciComp/Wavesurfer)
    %must be installed for this to work.
    
    clear vars
    
    d = dir('*.h5');
    str = {d.name};
    [s,~] = listdlg('PromptString','Select a file:', ...
                    'ListSize', [300, 600], ...
                    'SelectionMode','single', ...
                    'ListString',str);
                
    fileName = char(str(s));
    
    %Loads the dialogue selected data file as a structure to be passed
    %into different functions and loads the sample rate from the header
    wsDataStruct = ws.loadDataFile(fileName);
    sampleRate = wsDataStruct.header.Acquisition.SampleRate;
    
    %Detects start sample of TTL pulse generated by command to LED
    %[startStimulus] = detectTTL(wsDataStruct);
    startStimulus = 10000;
    
    %Generates matrix for ch1 and ch2 of raw data and plots it relative to
    %the standard deviation and mean, useful for checking if filter was
    %overzealous
    [ch1RawData,ch2RawData] = buildDataMatrix(wsDataStruct);
    [ch1Mean,ch2Mean,ch1StDev,ch2StDev] = getStats(ch1RawData,ch2RawData);
    plotData(ch1RawData,ch2RawData, ch1Mean,ch2Mean,ch1StDev,ch2StDev, sampleRate);
    
    %Filters out spikes and other abnormal activity in window that begins
    %with TTL pulse generating LED signal and ends after signal returns to
    %baseline [WIP, currently just plots until 500 mS after 0.5 seconds
    %into sweep
    ipscMaxDuration = 90;
    [filteredDataCh1,filteredDataCh2] = filterData(ch1RawData, ...
            ch2RawData, ch1Mean, ch2Mean, ch1StDev, ch2StDev, sampleRate, ...
            startStimulus, ipscMaxDuration);
        
    [ch1FilteredMean,ch2FilteredMean,ch1FilteredStDev,ch2FilteredStDev] = ...
        getStats(filteredDataCh1,filteredDataCh2);
    
    plotData(filteredDataCh1,filteredDataCh2, ch1FilteredMean, ...
        ch2FilteredMean,ch1FilteredStDev,ch2FilteredStDev, sampleRate);
    
    [restCurrent, endIPSC] = detectBaseline(filteredDataCh1, ...
        filteredDataCh2, ch1FilteredMean, ch2FilteredMean,...
        ch1FilteredStDev,ch2FilteredStDev, sampleRate, ...
        startStimulus);
        
    
end

function [ch1DataMatrix, ch2DataMatrix] = buildDataMatrix(inputData)
    %This function loads the channels (1 and 2) from the data structure
    %output from ws.loadDataFile. It first preallocates the matrix by
    %ensuring the last run was completed, then returns two matrices, one
    %for each channel. Since channel 1 and 2 are acquired simultaenously,
    %there is no need to test the length of channel 2. Channel 1 =
    %Headstage 1, Channel 2 = Headstage 2
    
    fields = fieldnames(inputData);
    
    firstSweepLength = size(inputData.(fields{2}).analogScans(:,1),1);
    lastSweepLength = size(inputData.(fields{end}).analogScans(:,1),1);
    
    if firstSweepLength == lastSweepLength;
        ch1DataMatrix = zeros(firstSweepLength, size(fields,1)-1);
        ch2DataMatrix = zeros(firstSweepLength, size(fields,1)-1);
    else
        ch1DataMatrix = zeros(firstSweepLength, size(fields,1)-2);
        ch2DataMatrix = zeros(firstSweepLength, size(fields,1)-2);
    end

    for fieldIdx = 2:(size(ch1DataMatrix,2)+1);
        ch1DataMatrix(:,fieldIdx-1) = inputData.(fields{fieldIdx}).analogScans(:,1);
        ch2DataMatrix(:,fieldIdx-1) = inputData.(fields{fieldIdx}).analogScans(:,2);
    end
end

function [ch1MeanVector,ch2MeanVector,ch1StdDevVector,ch2StdDevVector] = ...
                                    getStats(ch1DataMatrix,ch2DataMatrix)
    %This computes the mean and standard deviation across the rows of each
    %matrix, returning the average of all recordings and the standard
    %deviation for each time point
    
    ch1MeanVector = mean(ch1DataMatrix,2);
    ch1StdDevVector = std(ch1DataMatrix,0,2);
    
    ch2MeanVector = mean(ch2DataMatrix,2);
    ch2StdDevVector = std(ch2DataMatrix,0,2);
end
    
function [] = plotData(ch1DataMatrix,ch2DataMatrix, ch1MeanVector, ...
              ch2MeanVector, ch1StdDevVector, ch2StdDevVector, sampleRate)
    figure()
    xs =(1:size(ch1MeanVector))'/sampleRate;
    
    subplot(2,2,1)
    plotAverage(ch1MeanVector, ch1StdDevVector, sampleRate)
    
    subplot(2,2,2)
    plot(xs(1:5:end), ch1DataMatrix(1:5:end,:))
    
    subplot(2,2,3)
    plotAverage(ch2MeanVector, ch2StdDevVector, sampleRate)
    
    subplot(2,2,4)
    plot(xs(1:5:end), ch2DataMatrix(1:5:end,:))
    
    function [] = plotAverage(meanVector, stDevVector, sampleRate)
        x =(1:size(meanVector))'/sampleRate;
        lo = meanVector - stDevVector;
        hi = meanVector + stDevVector;
        hp = patch([x; x(end:-1:1); x(1)], [lo; hi(end:-1:1); lo(1)], 'r');
        
        hold on;
        hl = line(x,meanVector);
        set(hp, 'facecolor', [1 0.8 0.8], 'edgecolor', 'none');
        set(hl, 'color', 'r');
    end
    
end

function [filteredCh1,filteredCh2] = filterData(ch1DataMatrix, ...
    ch2DataMatrix, ch1MeanVector,ch2MeanVector,ch1StdDevVector, ...
    ch2StdDevVector, sampleRate, startStim, returnBaselineMS)

    endStim = startStim + returnBaselineMS*sampleRate/1000;

    [~,p] = size(ch1DataMatrix);
    ch1MeanMatrix = repmat(ch1MeanVector(startStim:endStim, :), 1, p);
    ch2MeanMatrix = repmat(ch2MeanVector(startStim:endStim, :), 1, p);
    
    ch1SigmaMatrix = repmat(ch1StdDevVector(startStim:endStim,:), 1, p);
    ch2SigmaMatrix = repmat(ch2StdDevVector(startStim:endStim,:), 1, p);
    
    ch1Outliers = abs(ch1DataMatrix(startStim:endStim, :) - ch1MeanMatrix) > 3*ch1SigmaMatrix;
    ch2Outliers = abs(ch2DataMatrix(startStim:endStim, :) - ch2MeanMatrix) > 3*ch2SigmaMatrix;

    filterCh1 = any(ch1Outliers,1);
    filterCh2 = any(ch2Outliers,1);
    commonFilter = logical(filterCh1 + filterCh2);

    filteredCh1 = ch1DataMatrix(:, logical(1-commonFilter));
    filteredCh2 = ch2DataMatrix(:, logical(1-commonFilter));
    
    %[~,I] = min(responseWindowCh1);
    %B = sortrows(responseWindowCh1',I)';
    %h = surface(B);
    %set(h,'LineStyle', 'none')
end

function [startTTL] = detectTTL(inputData)
    thresh = 2.5;
    
    fields = fieldnames(inputData);
    
    firstSweepLength = size(inputData.(fields{2}).analogScans(:,1),1);
    lastSweepLength = size(inputData.(fields{end}).analogScans(:,1),1);
    
    if firstSweepLength == lastSweepLength;
        TTLMatrix = zeros(firstSweepLength, size(fields,1)-1);
    else
        TTLMatrix = zeros(firstSweepLength, size(fields,1)-2);
    end

    for fieldIdx = 2:(size(TTLMatrix,2)+1);
        TTLMatrix(:,fieldIdx-1) = inputData.(fields{fieldIdx}).analogScans(:,6);
    end
    
    startTTL = find(TTLMatrix > thresh, 1);

end

function [baseline, endIPSC] = detectBaseline (ch1DataMatrix, ... 
        ch2DataMatrix, ch1MeanVector, ch2MeanVector, ch1StdDevVector, ...
        ch2StdDevVector, sampleRate, startStim)
    
    baseline = -50;
    endIPSC = 12000;
    
    assignin('base', 'ch1DataMatrix', ch1DataMatrix)
    assignin('base', 'ch2DataMatrix', ch2DataMatrix)
    assignin('base', 'ch1MeanVector', ch1MeanVector)
    assignin('base', 'ch2MeanVector', ch2MeanVector)
    assignin('base', 'ch2MeanVector', ch2MeanVector)
    
    assignin('base', 'ch1StdDevVector', ch1StdDevVector)
    assignin('base', 'ch2StdDevVector', ch2StdDevVector)
    
    assignin('base', 'startStim', startStim)
    assignin('base', 'baseline', baseline)
    assignin('base', 'endIPSC', baseline)
        
end