Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

MEAPatchAnalysis/MEAPatchAnalysis_010315.bak

Go to file
 
 
Cannot retrieve contributors at this time
1347 lines (1131 sloc) 74.5 KB
Raw Blame
Open in GitHub Desktop
classdef MEAPatchAnalysis < handle
properties
ParFilename=[];
plotParams=[];
currentDir=[];
rootDataDir=[];
overwrite;
sendEMailMessages=false;
email4Messages='shein.mark@gmail.com';
end
properties(SetAccess = protected)
Fs_ms
resamplingFactor
Par=[];
end
properties (Constant)
serverRootWin='\\storage.laur.corp.brain.mpg.de';
serverRootLinux='/storage/laur';
end
methods
function obj=MEAPatchAnalysis(rootDataDir) %class constructor
if nargin==0
if ispc
obj.rootDataDir=[obj.serverRootWin filesep' 'Data' filesep 'HembergerMike' filesep 'MEA'];
elseif isunix
obj.rootDataDir=[obj.serverRootLinux filesep' 'Data' filesep 'HembergerMike' filesep 'MEA'];
end
disp(['No root directory entered, default set to Root directory: ' num2str(obj.rootDataDir)]);
else
obj.rootDataDir=rootDataDir;
end
obj.currentDir=cd;
if ~strcmp(obj.currentDir,obj.rootDataDir)
disp(['Notice that you are not working in the default directory (this may change processing) : ' obj.rootDataDir]);
end
%initialize overwrite structure
obj.overwrite.spikeSorting=false;
obj.overwrite.patchExtraction=false;
obj.overwrite.STAExtraction=false;
obj.overwrite.patchTriggeredMEAExtraction=false;
end
function obj=getExcelData(obj,excelRecordingDataFileName)
Par.all.fieldsToExtract={'Exclude','Turtle','Date','Neuron','MEA file','Patch file','Series','Clamp','Vm','delay2Trig_ms','protocol','slab','mea','Sampling Correction','Heka Ch'};
Par.all.FsPatch=10e3; %for patch data
Par.all.FsMEA=20e3;
Par.all.patchDataUpsamplingFactor=Par.all.FsMEA/Par.all.FsPatch;
Par.all.patchSpikeTriggeredMEAData=1;
Par.all.prePatchSpikeOnMEA=20;%ms
Par.all.postPatchSpikeOnMEA=100;%ms
Par.all.preMEASpikeOnPatch=200;%ms
Par.all.postMEASpikeOnPatch=500;%ms
Par.all.patch2MEAtimeFac=0.997305;
Par.all.medFiltWin=15/1000*Par.all.FsMEA;
Par.all.medFiltDev=10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get data from excel file
defaultExcelName='MEA_DataSheet.xlsx';
if nargin==1
[num,txt,raw] = xlsread(defaultExcelName,1);
disp(['Extracting parameter from default excel sheet: ' defaultExcelName]);
else
[num,txt,raw] = xlsread(excelRecordingDataFileName,1);
end
%find valid table headers in excel data shit
pValidFields=logical(cell2mat(cellfun(@(x) sum(~isnan(x)),raw(1,:),'UniformOutput',0))); %When there is no title there will be a NaN value
recordingDataFieldNames=raw(1,pValidFields);
%remove excluded experiments
pFieldTmp=find(strcmp(recordingDataFieldNames,'Exclude'));
pExcludedExperiments=cell2mat(raw(2:end,pFieldTmp));
raw(find(pExcludedExperiments==1)+1,:)=[];
%find specific field data
for i=1:numel(Par.all.fieldsToExtract)
pFieldTmp=find(strcmp(recordingDataFieldNames,Par.all.fieldsToExtract{i}));
fieldShortName{i}=regexprep(Par.all.fieldsToExtract{i},'[^\w'']','');
Par.ExcelTable.( fieldShortName{i})=raw(2:end,pFieldTmp);
end
%Assemble data into separate experiment and sessions
pExperiments=logical(cell2mat(cellfun(@(x) sum(~isnan(x)),Par.ExcelTable.Turtle,'UniformOutput',0)));
startRawPerExperiment=find(pExperiments);
endRawPerExperiment=[startRawPerExperiment(2:end)-1;numel(pExperiments)];
allFields=fields(Par.ExcelTable);
for i=1:numel(startRawPerExperiment)
%define experiment name
Par.experimentNames{i}=[Par.ExcelTable.Turtle{startRawPerExperiment(i)} '_' num2str(Par.ExcelTable.Date{startRawPerExperiment(i)})];
%find the Neurons that were recorded in this specific experiment and check if they should not be excluded
tmpExperimentNeuron=Par.ExcelTable.Neuron(startRawPerExperiment(i):endRawPerExperiment(i));
tmpExperimentNeuronExclude=cell2mat(Par.ExcelTable.Exclude(startRawPerExperiment(i):endRawPerExperiment(i)));
%go over all Neurons in the experiment and callect all fields
NeuronNumberExp=cell2mat(tmpExperimentNeuron);
NeuronNumber=unique(NeuronNumberExp); %a session is defined by a patch in one Neuron which usually corresponds to not moving the grid possition
for j=1:numel(Par.all.fieldsToExtract)
tmpFieldValuesField=Par.ExcelTable.(allFields{j})(startRawPerExperiment(i):endRawPerExperiment(i));
for k=1:numel(NeuronNumber)
Par.(Par.experimentNames{i}).(['Neuron' num2str(NeuronNumber(k))]).(fieldShortName{j})=tmpFieldValuesField(find(NeuronNumber(k)==NeuronNumberExp));
end
end
end
%prepare other fields for the rest of the analysis
Exp=Par.experimentNames;
nExp=numel(Exp);
for i=1:nExp %go over all experiments
Neurons=fields(Par.(Exp{i}));
for j=1:numel(Neurons) %go over all Neurons within an experiment
tmpPatchDataName=cellfun(@(x) [obj.currentDir filesep Exp{i} filesep num2str(x) '.mat'],Par.(Exp{i}).(Neurons{j}).Patchfile,'UniformOutput',0);
Par.(Exp{i}).(Neurons{j}).fullPatchDataFilename=tmpPatchDataName;
MEATriggeredPatchFilename=[obj.currentDir filesep Exp{i} filesep Neurons{j} '_PatchSTA.mat'];
Par.(Exp{i}).(Neurons{j}).MEATriggeredPatchFilename=MEATriggeredPatchFilename;
seriesSortingName=[obj.currentDir filesep Exp{i} filesep Neurons{j} '_SpikeSortingData.mat'];
Par.(Exp{i}).(Neurons{j}).seriesSortingName=seriesSortingName;
spikeCrossCorrFileName=[obj.currentDir filesep Exp{i} filesep Neurons{j} '_spikeCrossCorr.mat'];
Par.(Exp{i}).(Neurons{j}).spikeCrossCorrFileName=spikeCrossCorrFileName;
patchSeries=Par.(Exp{i}).(Neurons{j}).Series;
nSeries=numel(patchSeries);
for k=1:nSeries %go over series in patch data
%check if multiple mcd files are recorded in the series
pSeperator4McdFiles=strfind(Par.(Exp{i}).(Neurons{j}).MEAfile{k},',');
if ~isempty(pSeperator4McdFiles)
Par.(Exp{i}).(Neurons{j}).MEAfile{k}=strsplit(Par.(Exp{i}).(Neurons{j}).MEAfile{k},',');
Par.(Exp{i}).(Neurons{j}).fullMEAFilenames{k}=cellfun(@(x) [obj.currentDir filesep Exp{i} filesep x '.mcd'],Par.(Exp{i}).(Neurons{j}).MEAfile{k},'UniformOutput',0);
else
Par.(Exp{i}).(Neurons{j}).fullMEAFilenames{k}={[obj.currentDir filesep Exp{i} filesep Par.(Exp{i}).(Neurons{j}).MEAfile{k} '.mcd']};
end
patchUpsampledFilename=[obj.currentDir filesep Exp{i} filesep Neurons{j} '_Ser' num2str(patchSeries{k}) '_UpsampledPatchData.mat'];
Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{k}=patchUpsampledFilename;
patchTriggeredMEAFilename=[obj.currentDir filesep Exp{i} filesep Neurons{j} '_Ser' num2str(patchSeries{k}) '_PatchTriggeredMEAData.mat'];
Par.(Exp{i}).(Neurons{j}).patchTriggeredMEAFilename{k}=patchTriggeredMEAFilename;
end
shortMEAfileNames=Par.(Exp{i}).(Neurons{j}).MEAfile{1};
if iscell(shortMEAfileNames) %separate for the case a series has more than one mcd file
Par.(Exp{i}).(Neurons{j}).sortingFileName=[obj.currentDir filesep Exp{i} filesep shortMEAfileNames{1} '_spikeSort' filesep 'SS_' shortMEAfileNames{1}];
else
Par.(Exp{i}).(Neurons{j}).sortingFileName=[obj.currentDir filesep Exp{i} filesep shortMEAfileNames '_spikeSort' filesep 'SS_' shortMEAfileNames];
end
end
end
obj.Par=Par;
end %getExcelData
function obj=plotPatchTriggeredMEAPhyiscalSpace(obj,h,exp,neuron,series,preSpike_ms,postSpike_ms)
Fs=obj.Par.all.FsMEA;
preSpike=preSpike_ms/1000*Fs;
postSpike=postSpike_ms/1000*Fs;
preSpikeData=obj.Par.all.prePatchSpikeOnMEA/1000*Fs;
pInData=(preSpikeData-preSpike+1):(preSpikeData+postSpike);
pSeries=find(cell2mat(obj.Par.(exp).(neuron).Series)==series);
if isempty(h)
f=figure('position',[100 100 800 800]);
h=axes;
end
if exist(obj.Par.(exp).(neuron).patchTriggeredMEAFilename{pSeries},'file')
load(obj.Par.(exp).(neuron).patchTriggeredMEAFilename{pSeries});
if isempty(MAll)
disp(['Cant plot - activity matrix in ' obj.Par.(exp).(neuron).patchTriggeredMEAFilename{pSeries} 'is empty']);
delete(h);
return;
end
else
disp(['Cant plot - non existing file: ' obj.Par.(exp).(neuron).patchTriggeredMEAFilename{pSeries}]);
delete(h);
return;
end
load(['layout_' obj.Par.(exp).(neuron).mea{pSeries}]);
%removes these lines after rerunnnig analysis
%channelNumbers=[1:14 16:60];
%MAll=MAll(:,:,pInData);
%MAll=bsxfun(@minus,MAll,mean(mean(MAll,2),3));
[hPlot, scale]=activityTracePhysicalSpacePlot(h,channelNumbers',squeeze(mean(MAll(1:end-1,:,:),2)),En,'DrawElectrodeNumbers',0);
[yE,xE]=size(En);
[hScaleBar]=addScaleBar(h,'xl_real',xE*[0 preSpike_ms+postSpike_ms],'yl_real',[0 yE*scale(2)]);
title([num2str(size(MAll,2)) ' spikes']);
end %plotPatchTriggeredMEAPhyiscalSpace
function obj=plotMEATriggeredPatchPhyiscalSpace(obj,h,exp,neuron,series,preSpike_ms,postSpike_ms)
Fs=obj.Par.all.FsMEA;
minNumberOfSpikes2Display=20;
DrawElectrodeNumbers=false;
preSpike=preSpike_ms/1000*Fs;
postSpike=postSpike_ms/1000*Fs;
preSpikeData=obj.Par.all.preMEASpikeOnPatch/1000*Fs;
pInData=(preSpikeData-preSpike+1):(preSpikeData+postSpike);
pSeries=find(cell2mat(obj.Par.(exp).(neuron).Series)==series);
if isempty(h)
f=figure('position',[100 100 800 800]);
h=axes;
end
if exist(obj.Par.(exp).(neuron).MEATriggeredPatchFilename,'file')
load(obj.Par.(exp).(neuron).MEATriggeredPatchFilename);
load(obj.Par.(exp).(neuron).seriesSortingName,'ic');
nSpikes=sum(nSpikesAll{pSeries},2);
else
disp(['Cant plot - non existing file: ' obj.Par.(exp).(neuron).patchTriggeredMEAFilename 'series ' num2str(series)]);
return;
end
load(['layout_' obj.Par.(exp).(neuron).mea{pSeries}]);
%removes these lines after rerunnnig analysis
pCh=(nSpikes>=minNumberOfSpikes2Display);
channelNumbers=ic{pSeries}(1,pCh)';
neuNumbers=ic{pSeries}(2,pCh);
STA=mMall{pSeries}(pCh,pInData);
STARand=mMallRand{pSeries}(pCh,pInData);
nSpikes=nSpikes(pCh);
STA=bsxfun(@minus,STA,mean(STA(:,1:preSpike),2));
%[hPlot,scale]=activityTracePhysicalSpacePlot(h,channelNumbers,STA,En,'scaling','noOverlap','DrawElectrodeNumbers',0,'LineWidth',1);
[hPlot,scale]=activityTracePhysicalSpacePlot(h,channelNumbers,STA,En,'DrawElectrodeNumbers',DrawElectrodeNumbers,'LineWidth',1,'scaling','noOverlap');
STARand=bsxfun(@minus,STARand,mean(STARand(:,1:preSpike),2));
[hPlot]=activityTracePhysicalSpacePlot(h,channelNumbers,STARand/scale,En,'scaling','none','traceColor',[0.8 0.2 0.2],'LineWidth',1,'DrawGrid',1);
%[hPlot]=activityTracePhysicalSpacePlot(h,channelNumbers,(STARand-scale(1))/scale(2),En,'scaling','none','traceColor',[0.8 0.2 0.2],'LineWidth',1,'DrawGrid',1);
[yE,xE]=size(En);
if strcmp(obj.Par.(exp).(neuron).Clamp{pSeries},'VC')
YUnitStr='pA';
elseif strcmp(obj.Par.(exp).(neuron).Clamp{pSeries},'CC')
YUnitStr='mV';
end
[hScaleBar]=addScaleBar(h,'xl_real',xE*[0 preSpike_ms+postSpike_ms],'yl_real',[0 yE*scale],'YUnitStr',YUnitStr);
%title([num2str(size(MAll,2)) ' spikes']);
end %plotMEATriggeredPatchPhyiscalSpace
function obj=plotMEATriggeredPatch(obj,f,exp,neuron,series,preSpike_ms,postSpike_ms)
Fs=obj.Par.all.FsMEA;
minNumberOfSpikes2Display=20;
preSpike=preSpike_ms/1000*Fs;
postSpike=postSpike_ms/1000*Fs;
preSpikeData=obj.Par.all.preMEASpikeOnPatch/1000*Fs;
pInData=(preSpikeData-preSpike+1):(preSpikeData+postSpike);
pSeries=find(cell2mat(obj.Par.(exp).(neuron).Series)==series);
if isempty(f)
f=figure('position',[100 100 1100 600]);
end
if exist(obj.Par.(exp).(neuron).MEATriggeredPatchFilename,'file')
load(obj.Par.(exp).(neuron).MEATriggeredPatchFilename);
load(obj.Par.(exp).(neuron).seriesSortingName,'ic');
nSpikes=sum(nSpikesAll{pSeries},2);
else
disp(['Cant plot - non existing file: ' obj.Par.(exp).(neuron).patchTriggeredMEAFilename 'series ' num2str(series)]);
return;
end
pCh=(nSpikes>=minNumberOfSpikes2Display);
channelNumbers=ic{pSeries}(1,pCh);
neuNumbers=ic{pSeries}(2,pCh);
STA=mMall{pSeries}(pCh,pInData);
STARand=mMallRand{pSeries}(pCh,pInData);
nSpikes=nSpikes(pCh);
if numel(STA)>0
STA=bsxfun(@minus,STA,mean(STA,2));
STARand=bsxfun(@minus,STARand,mean(STARand,2));
if strcmp(obj.Par.(exp).(neuron).Clamp{pSeries},'VC')
YUnitStr='I [pA]';
elseif strcmp(obj.Par.(exp).(neuron).Clamp{pSeries},'CC')
YUnitStr='V [mV]';
end
nSTAs=size(STA,1);
mMax=[3:15];
nMax=round(mMax*11/7);
fontSize=round(8:(-4/numel(mMax)):4);
pMN=find(mMax.*nMax>nSTAs,1,'first');
mMax=mMax(pMN);
nMax=nMax(pMN);
fontSize=fontSize(pMN);
yl=[-nanstd(STA(:))*5 nanstd(STA(:))*5];
for i=1:nSTAs
if pMN>12
h(i)=subaxis(mMax,nMax,i,'M',0.001,'SH',0.001,'SV',0.001);
else
h(i)=subaxis(mMax,nMax,i,'M',0.05,'SH',0.02,'SV',0.03);
end
plot(((1-preSpike):postSpike)*1000/Fs,STA(i,:));hold on;
plot(((1-preSpike):postSpike)*1000/Fs,STARand(i,:),'r');
xlim([-preSpike_ms postSpike_ms]);
ylim(yl);
line([0 0],yl,'color',[0.8 0.8 0.8]);
title(['Ch=' num2str(channelNumbers(i)) '-' num2str(neuNumbers(i)) ',N=' num2str(nSpikes(i))],'fontSize',fontSize);
end
set(h,'fontSize',fontSize);
set(h,'XTickLabel',[],'YTickLabel',[]);
xEdgePlots=max(nSTAs-nMax+1,1):nSTAs;
if numel(xEdgePlots)~=1
set(cell2mat(get(h(xEdgePlots),'XLabel')),'String','T [ms]'); %this does not support a single axis where get returns a double and not a cell array
set(h(xEdgePlots),'XTickLabelMode','auto');
else
xlabel(h(xEdgePlots),'T [ms]');
end
yEdgePlots=1:nMax:nSTAs;
if numel(yEdgePlots)~=1
set(cell2mat(get(h(yEdgePlots),'YLabel')),'String',YUnitStr); %this does not support a single axis where get returns a double and not a cell array
else
ylabel(h(yEdgePlots),YUnitStr);
end
set(h(yEdgePlots),'YTickLabelMode','auto');
end
end %plotMEATriggeredPatchPhyiscalSpace
function obj=runSpikeSorting(obj,Exp,Neurons,channels)
if nargin<4
channels=[];
end
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
fastPrinting=0;
if obj.overwrite.spikeSorting
overwriteSpikeExtraction=1;
overwriteFeatureExtraction=1;
overwriteClustering=1;
overwriteMerging=1;
overwriteFitting=1;
overwritePostProcessignAnalysis=1;
else
overwriteSpikeExtraction=0;
overwriteFeatureExtraction=0;
overwriteClustering=0;
overwriteMerging=0;
overwriteFitting=0;
overwritePostProcessignAnalysis=0;
end
try
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp
disp(['******Starting spike sorting analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
for j=1:numel(Neurons)
disp(['* ' num2str(Neurons{j})]);
c=1;
shortMEAfileNames={};
fullMEAfileNames={};
for m=1:numel(obj.Par.(Exp{i}).(Neurons{j}).MEAfile)
if iscell(obj.Par.(Exp{i}).(Neurons{j}).MEAfile{m}) %separate for the case a series has more than one mcd file
shortMEAfileNames(c:(c+numel(obj.Par.(Exp{i}).(Neurons{j}).MEAfile{m})-1))=obj.Par.(Exp{i}).(Neurons{j}).MEAfile{m};
fullMEAfileNames(c:(c+numel(obj.Par.(Exp{i}).(Neurons{j}).MEAfile{m})-1))=obj.Par.(Exp{i}).(Neurons{j}).fullMEAFilenames{m};
else
shortMEAfileNames{c}=obj.Par.(Exp{i}).(Neurons{j}).MEAfile{m};
fullMEAfileNames{c}=obj.Par.(Exp{i}).(Neurons{j}).fullMEAFilenames{m}{1};
end
c=numel(shortMEAfileNames)+1;
end
tmpSortingFileName=obj.Par.(Exp{i}).(Neurons{j}).sortingFileName;
if obj.overwrite.spikeSorting || ~exist([tmpSortingFileName '.mat'],'file')
if ispc
MCR=MCRackRecording(fullMEAfileNames');
elseif isunix
MCR=MCRackRecordingNeuroshare(fullMEAfileNames');
end
channelNumbers=MCR.channelNumbers;
if mean(strcmp(obj.Par.(Exp{i}).(Neurons{j}).mea,'40_Hexa'))
gridSize=9;
gridExt=3;
%The lines below should be removed once the layout is determined automatically by the MCRackRecording object
load('layout_40_Hexa');
MCR.chLayoutNumbers=En;
MCR.chLayoutNames=Ena;
%remove from layout channels not in the channel list
unrecordedChannels=setdiff(En(~isnan(En)),MCR.channelNumbers);
for k=1:numel(unrecordedChannels)
p=find(En==unrecordedChannels(k));
MCR.chLayoutNumbers(p)=NaN;
MCR.chLayoutNames{p}=[];
end
else
gridSize=3;
gridExt=1;
end
[sortingDir,name,ext] = fileparts(tmpSortingFileName);
if isempty(channels)
fileNames=spikeSortingMultiChannel(MCR,'sortingDir',sortingDir,'minSpikesTotal',2,'minSpikesCluster',1,'gridSize',gridSize,'gridExt',gridExt,'fastPrinting',fastPrinting,...
'overwriteSpikeExtraction',overwriteSpikeExtraction,'overwriteFeatureExtraction',overwriteFeatureExtraction,'overwriteClustering',overwriteClustering,'overwriteMerging',overwriteMerging,'overwriteFitting',overwriteFitting,'overwritePostProcessignAnalysis',overwritePostProcessignAnalysis);
else
fileNames=spikeSortingMultiChannel(MCR,'selectedChannels',channelNumbers,'sortingDir',sortingDir,'minSpikesTotal',2,'minSpikesCluster',1,'gridSize',gridSize,'gridExt',gridExt,'fastPrinting',fastPrinting,...
'overwriteSpikeExtraction',overwriteSpikeExtraction,'overwriteFeatureExtraction',overwriteFeatureExtraction,'overwriteClustering',overwriteClustering,'overwriteMerging',overwriteMerging,'overwriteFitting',overwriteFitting,'overwritePostProcessignAnalysis',overwritePostProcessignAnalysis);
end
Triggers=MCR.getTrigger;
load(fileNames.fittingFile);
load(fileNames.mergedAvgWaveformFile,'avgWF','isNoise');
save(tmpSortingFileName,'t','ic','Triggers','channelNumbers','avgWF','isNoise','-v7.3');
end
end
end
catch errorMsg
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['An error occured while running runSpikeSorting on ' getenv('COMPUTERNAME') ' session ' num2str(i) '/' num2str(nExp)],errorMsg.getReport);
end
rethrow(errorMsg);
end
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['runSpikeSorting completed on ' getenv('COMPUTERNAME') ,' session success: ' num2str(nExp)]);
end
end %runSpikeSorting
function [obj,patchData]=runPatchDataExtraction(obj,Exp,Neurons,Series)
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
if nargout==2
save2file=0;
patchData=[];
else
save2file=1;
end
loadedPatchFile=''; %to not load the patch file if loaded for the previous neuron
obj.Fs_ms=obj.Par.all.FsMEA/1000;
try
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting patch-MEA analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
%check for sampling issues in the MEA system in this experiment and update resamplingFactor
correctSampling=obj.Par.(Exp{i}).(Neurons{1}).SamplingCorrection{1};
if correctSampling
obj.resamplingFactor=obj.Par.all.patch2MEAtimeFac;
else
obj.resamplingFactor=1;
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' num2str(Neurons{j})]);
% Extract patch data (including upsampling
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
MEAfileNames=obj.Par.(Exp{i}).(Neurons{j}).fullMEAFilenames;
clamp=obj.Par.(Exp{i}).(Neurons{j}).Clamp;
delay2Trig=obj.Par.(Exp{i}).(Neurons{j}).delay2Trig_ms;
tmpPatchDataName=obj.Par.(Exp{i}).(Neurons{j}).fullPatchDataFilename;
%check if to only analyze specific series, else analyzes all series with exluded series where patch data is not available
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
pSeries(isnan(cell2mat(patchSeries(pSeries))))=[];
nSeries=numel(pSeries);
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
for k=1:nSeries %go over series in patch data
ind=pSeries(k);
disp(['Analyzing Series ' num2str(patchSeries{ind})]);
patchUpsampledFilename=obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{ind};
hekaCh=obj.Par.(Exp{i}).(Neurons{j}).HekaCh{ind};
if obj.overwrite.patchExtraction || ~exist(patchUpsampledFilename,'file')
%load patch data only if it was not loaded for the previous Neuron
if ~strcmp(loadedPatchFile,tmpPatchDataName{k})
load(tmpPatchDataName{k},'data'); %it is assumned that the same Neuron will always belong to the same patch recording
loadedPatchFile=tmpPatchDataName{k};
end
[ch , t]= parseHeka(data,patchSeries{ind},hekaCh);
units=ch.units;
nPatchSamples=numel(t.ms);
nPatchSweeps=size(ch.vm,2);
tMp=(1000/obj.Par.all.FsMEA/obj.resamplingFactor):(1000/obj.Par.all.FsMEA/obj.resamplingFactor):max(t.ms);
resampledSamples=numel(tMp);
Mp=zeros(nPatchSweeps,resampledSamples); %initialize data array
for m=1:nPatchSweeps %upsample data
Mp(m,:) = interp1(t.ms,ch.vm(:,m),tMp,'spline');
end
if save2file
save(patchUpsampledFilename,'Mp','tMp','nPatchSweeps','nPatchSamples','resampledSamples','units','-v7.3','-mat');
else
patchData.Mp=Mp;
patchData.tMp=tMp;
patchData.nPatchSweeps=nPatchSweeps;
patchData.nPatchSamples=nPatchSamples;
patchData.resampledSamples=resampledSamples;
patchData.units=units;
end
end
end %for loop over series
end %for loop over Neurons
end %for loop over recordings
catch errorMsg
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['An error occured while running runPatchDataExtractionAndMEASTA on ' getenv('COMPUTERNAME') ' session ' num2str(i) '/' num2str(nExp)],errorMsg.getReport);
end
rethrow(errorMsg);
end
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['runPatchDataExtractionAndMEASTA completed on ' getenv('COMPUTERNAME') ,' session success: ' num2str(nExp)]);
end
end %runPatchDataExtraction
function [obj,MEASTA]=runMEASTAs(obj,Exp,Neurons,Series)
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
spikePeakDetectionWinMs=3;
minFilterCrossMs=0.2;
minSpikePeakValue=0; %usually mV
if nargout==2
save2file=0;
patchData=[];
obj.overwrite.patchTriggeredMEAExtraction=1;
else
save2file=1;
end
loadedPatchFile=''; %to not load the patch file if loaded for the previous neuron
obj.Fs_ms=obj.Par.all.FsMEA/1000;
spikePeakDetectionWinSamples=spikePeakDetectionWinMs*obj.Fs_ms;
minFilterCrossSamples=minFilterCrossMs*obj.Fs_ms;
try
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting patch-MEA analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
%check for sampling issues in the MEA system in this experiment and update resamplingFactor
correctSampling=obj.Par.(Exp{i}).(Neurons{1}).SamplingCorrection{1};
if correctSampling
obj.resamplingFactor=obj.Par.all.patch2MEAtimeFac;
else
obj.resamplingFactor=1;
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' num2str(Neurons{j})]);
% Extract patch data (including upsampling
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
MEAfileNames=obj.Par.(Exp{i}).(Neurons{j}).fullMEAFilenames;
clamp=obj.Par.(Exp{i}).(Neurons{j}).Clamp;
delay2Trig=obj.Par.(Exp{i}).(Neurons{j}).delay2Trig_ms;
tmpPatchDataName=obj.Par.(Exp{i}).(Neurons{j}).fullPatchDataFilename;
%check if to only analyze specific series, else analyzes all series
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
pSeries(isnan(cell2mat(patchSeries(pSeries))))=[];
nSeries=numel(pSeries);
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
for k=1:nSeries %go over series in patch data
ind=pSeries(k);
disp(['Analyzing Series ' num2str(patchSeries{ind})]);
patchUpsampledFilename=obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{ind};
%extract patch triggered MEA traces
patchTriggeredMEAFilename=obj.Par.(Exp{i}).(Neurons{j}).patchTriggeredMEAFilename{ind};
if (strcmp(clamp{ind},'CC') && (delay2Trig{ind}==2000 || delay2Trig{ind}==3000)) && (~exist(patchTriggeredMEAFilename,'file') || obj.overwrite.patchTriggeredMEAExtraction) %calculate patch triggered MEA recording
%load patch data only if it was not loaded for the previous Neuron
if ~strcmp(loadedPatchFile,tmpPatchDataName{k})
load(tmpPatchDataName{k},'data'); %it is assumned that the same Neuron will always belong to the same patch recording
loadedPatchFile=tmpPatchDataName{k};
end
load(patchUpsampledFilename);
nSamples=(obj.Par.all.prePatchSpikeOnMEA+obj.Par.all.postPatchSpikeOnMEA)*obj.Fs_ms;
if ispc
MCR=MCRackRecording(MEAfileNames{ind});
elseif isunix
MCR=MCRackRecordingNeuroshare(MEAfileNames{ind});
end
Triggers=MCR.getTrigger;
channelNumbers=MCR.channelNumbers;
cumsumNSpikes=0;
nSpikes=[];
pSpikes={};
startTimesOnMcd={};
MAll=[];
%check for the case when Mike did not turn MEA recording on time and only the first trigger is missing
if nPatchSweeps~=numel(Triggers{1})
disp('Notice: number of sweeps and trigers is not identical!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
badRecording=true;
Triggers{1}=[NaN Triggers{1}];
else
badRecording=false;
end
for m=1:nPatchSweeps
if badRecording && m==1
pSpikes{1}=[];
startTimesOnMcd{1}=[];
nSpikes(1)=0;
else
med = fastmedfilt1d(Mp(m,:), obj.Par.all.medFiltWin,-fliplr(Mp(m,1:(obj.Par.all.medFiltWin/2))),-fliplr(Mp(m,end+1-(obj.Par.all.medFiltWin/2):end)))';
absMed=abs(Mp(m,:)-med);
medDev = fastmedfilt1d(absMed,obj.Par.all.medFiltWin,-fliplr(absMed(1:(obj.Par.all.medFiltWin/2))),-fliplr(absMed(end+1-(obj.Par.all.medFiltWin/2):end)))'*1.4826;
patchSpikesLogical=Mp(m,:)>med+obj.Par.all.medFiltDev*medDev;
pUpCross=find((patchSpikesLogical(2:end)-patchSpikesLogical(1:end-1))==1);
pDnCross=find((patchSpikesLogical(2:end)-patchSpikesLogical(1:end-1))==-1);
%check that the duration of med filter crossing is not too short
if numel(pDnCross)==(numel(pUpCross))
pUpCross((pDnCross-pUpCross)<minFilterCrossSamples)=[];
elseif numel(pDnCross)==(numel(pUpCross)-1)
pDnCross=[pDnCross numel(patchSpikesLogical)];
pUpCross((pDnCross-pUpCross)<minFilterCrossSamples)=[];
else
%does not remove samples from pUpCross
end
spikeLocal=zeros(spikePeakDetectionWinSamples,numel(pUpCross));
pPeakMp=bsxfun(@plus,pUpCross,(1:spikePeakDetectionWinSamples)');
pPeakMp(pPeakMp>resampledSamples)=resampledSamples;
spikeLocal(:)=Mp(m,pPeakMp);
[vMax,pMax]=max(spikeLocal);
pSpikes{m}=pUpCross+pMax;
pSpikes{m}(pSpikes{m}<minSpikePeakValue)=[];
nSpikes(m)=numel(pSpikes{m});
% A criterion with spike durations of at least 1ms should be added to false positive spikes
%plot(Mp(m,:));hold on;plot(patchSpikesLogical*3,'r');plot(med+obj.Par.all.medFiltDev*medDev,'g');plot(pSpikes{m},Mp(m,pSpikes{m}),'ok');
nCh=numel(MCR.channelNumbers);
M=zeros(nCh+1,numel(pSpikes{m}),nSamples);
patchTmp=zeros(nSpikes(m),nSamples);
p=bsxfun(@plus,pSpikes{m}',-(obj.Par.all.prePatchSpikeOnMEA*obj.Fs_ms):(obj.Par.all.postPatchSpikeOnMEA*obj.Fs_ms-1));
p(p<=0)=1;p(p>nPatchSamples)=nPatchSamples;
patchTmp(:)=Mp(m,p);
M(nCh+1,:,:)=patchTmp;
% ********This part should be update to be compatible with a general electrode array
%startTimes=round(20*(Triggers{1}(m)-delay2Trig{ind}* 0.997275+pSpikes{m}/obj.Fs_ms* 0.997275-obj.Par.all.prePatchSpikeOnMEA))/20;
startTimes=round(obj.Fs_ms*(Triggers{1}(m)-delay2Trig{ind}*obj.resamplingFactor+pSpikes{m}/obj.Fs_ms-obj.Par.all.prePatchSpikeOnMEA))/obj.Fs_ms;
M(1:nCh,:,:)=MCR.getData(MCR.channelNumbers,startTimes,obj.Par.all.prePatchSpikeOnMEA+obj.Par.all.postPatchSpikeOnMEA);
MAll(:,cumsumNSpikes+1:cumsumNSpikes+nSpikes(m),:)=M;
cumsumNSpikes=cumsumNSpikes+nSpikes(m);
spikeTimesOnMcd{m}=startTimes+obj.Par.all.prePatchSpikeOnMEA;
end
end
if ~isempty(MAll)
MAll(1:nCh,:,:)=bsxfun(@minus,MAll(1:nCh,:,:),mean(mean(MAll(1:nCh,:,:),2),3)); %remove average
end
if save2file
save(patchTriggeredMEAFilename,'MAll','nSpikes','pSpikes','channelNumbers','spikeTimesOnMcd','-v7.3','-mat');
else
MEASTA.MAll=MAll;
MEASTA.nSpikes=nSpikes;
MEASTA.channelNumbers=channelNumbers;
MEASTA.spikeTimesOnMcd=spikeTimesOnMcd;
end
end
end %for loop over series
end %for loop over Neurons
end %for loop over recordings
catch errorMsg
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['An error occured while running runPatchDataExtractionAndMEASTA on ' getenv('COMPUTERNAME') ' session ' num2str(i) '/' num2str(nExp)],errorMsg.getReport);
end
rethrow(errorMsg);
end
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['runPatchDataExtractionAndMEASTA completed on ' getenv('COMPUTERNAME') ,' session success: ' num2str(nExp)]);
end
end %runMEASTAs
function [patchTraces,tPatchTraces,patchUnits]=getResampledPatchTraces(obj,Exp,Neurons,Series)
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
loadedPatchFile='';
obj.Fs_ms=obj.Par.all.FsMEA/1000;
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting patch-MEA analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' num2str(Neurons{j})]);
% Extract patch data (including upsampling
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
clamp=obj.Par.(Exp{i}).(Neurons{j}).Clamp;
delay2Trig=obj.Par.(Exp{i}).(Neurons{j}).delay2Trig_ms;
%check if to only analyze specific series, else analyzes all series
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
pSeries(isnan(cell2mat(patchSeries(pSeries))))=[];
nSeries=numel(pSeries);
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
for k=1:nSeries %go over series in patch data
disp(['Analyzing Series ' num2str(pSeries(k))]);
patchUpsampledFilename=obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{pSeries(k)};
load(patchUpsampledFilename,'units','Mp','tMp');
tPatchTraces{i,j,k}=tMp*obj.resamplingFactor;
patchTraces{i,j,k}=Mp;
patchUnits{i,j,k}=units;
end %for loop over series
end %for loop over Neurons
end %for loop over recordings
end %getResampledPatchTraces
function [patchSpikes,MEAspikes]=getMEATimes4PatchSpikes(obj,MEANeuron,Exp,Neurons,Series)
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
loadedPatchFile='';
obj.Fs_ms=obj.Par.all.FsMEA/1000;
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting patch-MEA analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' num2str(Neurons{j})]);
% Extract patch data (including upsampling
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
clamp=obj.Par.(Exp{i}).(Neurons{j}).Clamp;
delay2Trig=obj.Par.(Exp{i}).(Neurons{j}).delay2Trig_ms;
%check if to only analyze specific series, else analyzes all series
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
pSeries(isnan(cell2mat(patchSeries(pSeries))))=[];
nSeries=numel(pSeries);
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
patchSpikes=cell(1,nSeries);
MEAspikes=cell(1,nSeries);
for k=1:nSeries %go over series in patch data
disp(['Analyzing Series ' num2str(pSeries(k))]);
patchUpsampledFilename=obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{pSeries(k)};
load(patchUpsampledFilename,'units','nPatchSweeps','nPatchSamples');
seriesSortingName=obj.Par.(Exp{i}).(Neurons{j}).seriesSortingName;
load(seriesSortingName,'Triggers','t','ic');
neu=find(ic{pSeries(k)}(1,:)==MEANeuron(1) & ic{pSeries(k)}(2,:)==MEANeuron(2));
tTmp=t{pSeries(k)}(ic{pSeries(k)}(3,neu):ic{pSeries(k)}(4,neu));
%icTmp=[1;1;1;numel(tTmp)];
%extract patch triggered MEA traces
patchTriggeredMEAFilename=obj.Par.(Exp{i}).(Neurons{j}).patchTriggeredMEAFilename{pSeries(k)};
if ( strcmp(clamp{pSeries(k)},'CC') && (delay2Trig{pSeries(k)}==2000 || delay2Trig{pSeries(k)}==3000) )
load(patchTriggeredMEAFilename,'pSpikes');
patchSpikes{k}=cell(1,nPatchSweeps);
MEAspikes{k}=cell(1,nPatchSweeps);
%check for the case when Mike did not turn MEA recording on time and only the first trigger is missing
for m=1:nPatchSweeps
patchSpikes{k}{m}=pSpikes{m}/obj.Fs_ms;
startTime=Triggers{pSeries(k)}(m);%+delay2Trig{pSeries(k)}*obj.resamplingFactor;
MEAspikes{k}{m}=tTmp(find( tTmp> startTime & tTmp<(startTime+nPatchSamples/obj.Fs_ms) ))-startTime;
%plot(MEAspikes{k}{m},ones(1,numel(MEAspikes{k}{m})),'.');hold on;plot(patchSpikes{k}{m},ones(1,numel(patchSpikes{k}{m})),'or');
end
%plot(tTmp,ones(1,numel(tTmp)),'.');hold on;T=(Triggers{pSeries(k)}+delay2Trig{pSeries(k)});plot(T,ones(1,numel(T)),'or');plot(T+nPatchSamples/obj.Fs_ms,ones(1,numel(T)),'og');
else
patchSpikes{k}=[];
MEAspikes{k}=[];
end
end %for loop over series
end %for loop over Neurons
end %for loop over recordings
end %getMEATimes4PatchSpikes
function [obj,CC]=runSpikeCrossCorr(obj,Exp,Neurons)
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
if nargout>1
save2file=false;
else
save2file=true;
end
win=1000*60*1;
sigma=3;
kern=fspecial('gaussian',[1 sigma*8],sigma);
lag=200;
obj.Fs_ms=obj.Par.all.FsMEA/1000;
try
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting spike cross corr analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' num2str(Neurons{j})]);
load(obj.Par.(Exp{i}).(Neurons{j}).sortingFileName);
nNeu=size(ic,2);
%{
startTimes=round(max(0,min(t)-1000)):win:round(max(t)+1000);
nStartTimes=numel(startTimes);
mX=zeros(nStartTimes,nNeu,nNeu);
pmX=zeros(nStartTimes,nNeu,nNeu);
for c=1:nStartTimes
M=squeeze(ConvBurstMatrix(BuildBurstMatrix(ic,round(t),startTimes(c),win),kern))';
[Xtmp,lags] = xcorr(M,lag,'coeff');
[mXtmp,pmXtmp] = nanmax(Xtmp,[],1);
mX(c,:,:)=reshape(mXtmp,[nNeu nNeu]);
pmX(c,:,:)=reshape(pmXtmp,[nNeu nNeu]);
end
%}
%use cconv to calculate cross-correlation
%initialize parameters
mX=zeros(nNeu,nNeu);
pmX=zeros(nNeu,nNeu);
X=zeros(nNeu,nNeu,lag*2+1);
%M=squeeze(ConvBurstMatrix(BuildBurstMatrix(ic,round(t),0,round(max(t))),kern));
for m=1:nNeu
M1=convn(squeeze(BuildBurstMatrix([1;1;1;ic(4,m)-ic(3,m)+1],round(t(ic(3,m):ic(4,m))),0,round(max(t)))),kern,'same');
for n=(m+1):nNeu
M2=convn(squeeze(BuildBurstMatrix([1;1;1;ic(4,n)-ic(3,n)+1],round(t(ic(3,n):ic(4,n))),0,round(max(t)))),kern,'same');
[X(m,n,:)] = xcorr(M1,M2,lag,'coeff');
[mX(m,n),pmX(m,n)] = nanmax(X(m,n,:));
end
end
if save2file
save(obj.Par.(Exp{i}).(Neurons{j}).spikeCrossCorrFileName,'X','mX','pmX');
else
CC.X{i,j}=X;
CC.mX{i,j}=mX;
CC.pmX{i,j}=pmX;
end
%{
load(obj.Par.(Exp{i}).(Neurons{j}).seriesSortingName);
nSeries=numel(patchSeries);
for k=1:nSeries %go over series in patch data
disp(['Analyzing Series ' num2str(k)]);
load(obj.Par.(Exp{i}).(Neurons{j}).patchTriggeredMEAFilename{k},'pSpikes');
end %for loop over series
%}
end %for loop over Neurons
end %for loop over recordings
catch errorMsg
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['An error occured while running runSpikeCrossCorrAnalysis on ' getenv('COMPUTERNAME') ' session ' num2str(i) '/' num2str(nExp)],errorMsg.getReport);
end
rethrow(errorMsg);
end
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['runPatchDataExtractionAndMEASTA completed on ' getenv('COMPUTERNAME') ,' session success: ' num2str(nExp)]);
end
end %runPatchDataExtractionAndAnalysis
function [obj,triggeredData]=runMEATriggeredAvgOnPatchData(obj,Exp,Neurons,Series,MEAneuron)
triggeredData=[];
if isempty(obj.ParFilename)
error('Can not execute method without parameters property (''Par'')');
end
if nargout==2
calculateFullTraces=true;
else
calculateFullTraces=false;
end
obj.Fs_ms=obj.Par.all.FsMEA/1000;
%calculating spike triggered averages
preSpikeSamples=obj.Par.all.preMEASpikeOnPatch*obj.Fs_ms;
postSpikeSamples=obj.Par.all.postMEASpikeOnPatch*obj.Fs_ms;
winSpikeSamples=preSpikeSamples+postSpikeSamples;
try
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Starting patch-MEA analysis for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
%check for sampling issues in the MEA system in this experiment and update resamplingFactor
correctSampling=obj.Par.(Exp{i}).(Neurons{1}).SamplingCorrection{1};
if correctSampling
obj.resamplingFactor=obj.Par.all.patch2MEAtimeFac;
else
obj.resamplingFactor=1;
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' Neurons{j}]);
if calculateFullTraces || obj.overwrite.STAExtraction || ~exist(obj.Par.(Exp{i}).(Neurons{j}).seriesSortingName,'file') || ~exist(obj.Par.(Exp{i}).(Neurons{j}).MEATriggeredPatchFilename,'file')
load(obj.Par.(Exp{i}).(Neurons{j}).sortingFileName); %load spike sorting data
if nargin<5
tAll=t;
icAll=ic;
TriggersAll=Triggers;
else
if size(MEAneuron{1},2)==1 %change to check if all channels are valid
[tAll,icAll]=RemainNeurons(t,ic,MEAneuron{:});
TriggersAll=Triggers;
else
error('Invalid MEA neuron selected!!!');
end
end
%check if to only analyze specific series, else analyzes all series
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
pSeries(isnan(cell2mat(patchSeries(pSeries))))=[];
nSeries=numel(pSeries);
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
%examine all patch series and get triggers per series
nSweepsAll=zeros(1,nSeries);
for s=1:nSeries
tmp=load(obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{pSeries(s)},'nPatchSweeps'); %load upsampled patch data
nSweepsAll(s)=tmp.nPatchSweeps;
end
endTrigSweep=cumsum(nSweepsAll);
startTrigSweep=[1 endTrigSweep(1:end-1)+1];
if numel(TriggersAll{1})~=endTrigSweep(end)
warning(['The number of triggers (' num2str(numel(TriggersAll{1})) ') is different from the number of patch sweeps(' num2str(endTrigSweep(end)) ')!!!!! Check that mcd files are not missing!']);
%check for the case when Mike did not turn MEA recording on time and only the first trigger is missing
if numel(TriggersAll{2})==endTrigSweep(end)
TrigOnOffDelay=[TriggersAll{2}(1:end-1)-TriggersAll{1} TriggersAll{2}(2:end)-TriggersAll{1}];
minTrigOnOffDelay=min(TrigOnOffDelay(find(TrigOnOffDelay>0)));
TriggersAll{1}=TriggersAll{2}-minTrigOnOffDelay;
disp(['Missing trigger onsets were replaced by using a deviation of ' num2str(minTrigOnOffDelay) 'ms from the trigger offset'])
else
disp([ Exp{i} '-' Neurons{j} ' not analyzed due to trigger missmatch']);
continue;
end
end
delay2Trig=obj.Par.(Exp{i}).(Neurons{j}).delay2Trig_ms;
mMall=cell(1,nSeries);
mMallRand=cell(1,nSeries);
nSpikesAll=cell(1,nSeries);
nSpikesRandAll=cell(1,nSeries);
Triggers=cell(1,nSeries);
t=cell(1,nSeries);
tS=cell(1,nSeries);
ic=cell(1,nSeries);
tStart=zeros(1,nSeries);
for k=1:nSeries %go over series in patch data
ind=pSeries(k);
%separate spike sorting according to patch series
load(obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{ind},'Mp','nPatchSweeps'); %load upsampled patch data
%{
filterData=1;
if filterData
%[z1,p1,k1] = butter(8,150/20e3,'low');
[z1,p1,k1] = butter(8,[190 210]/20e3,'stop');
[sos1,G1] = zp2sos(z1,p1,k1);
Mp=filtfilt(sos1,G1,Mp')';
[z1,p1,k1] = butter(8,[390 410]/20e3,'stop');
[sos1,G1] = zp2sos(z1,p1,k1);
Mp=filtfilt(sos1,G1,Mp')';
end
%}
Triggers{k}=TriggersAll{1}(startTrigSweep(ind):endTrigSweep(ind))-delay2Trig{ind}*obj.resamplingFactor;%get triggers and
%Triggers marks the beginning of the patch recording
nPatchSamples=size(Mp,2);
[t{k},ic{k},tStart(k)]=CutSortChannel(tAll,icAll,Triggers{k}(1),Triggers{k}(end)+nPatchSamples/obj.Fs_ms,false);
tS{k}=t{k}(randperm(numel(t{k})));
Triggers{k}=Triggers{k}-tStart(k);
nNeurons=size(ic{k},2);
mMall{k}=zeros(nNeurons,winSpikeSamples);
mMallRand{k}=zeros(nNeurons,winSpikeSamples);
nSpikesAll{k}=zeros(nNeurons,nPatchSweeps);
nSpikesRandAll{k}=zeros(nNeurons,nPatchSweeps);
for m=1:nNeurons
%disp(['starting Neuron ' num2str(m)]);
tTmp=t{k}(ic{k}(3,m):ic{k}(4,m));
icTmp=[1;1;1;numel(tTmp)];
tTmpS=tS{k}(ic{k}(3,m):ic{k}(4,m));
%plot(tTmp,ones(1,numel(tTmp)),'.');hold on;T=Triggers{k};plot(T,ones(1,numel(T)),'or');
%Msort=BuildBurstMatrix(icSeriesTmp,round(tSeriesTmp*obj.Fs_ms),round(({1}-delay2Trig{k})*obj.Fs_ms),nPatchSamples);
Msort=BuildBurstMatrix(icTmp,round(tTmp*obj.Fs_ms),round(Triggers{k}*obj.Fs_ms),nPatchSamples);
MsortS=BuildBurstMatrix(icTmp,round(tTmpS*obj.Fs_ms),round(Triggers{k}*obj.Fs_ms),nPatchSamples);
%MShuffled=zeros(size(Msort)); of shuffling
Mall=[];
MallRand=[];
spikesToRemoveAll=[];
spikesToRemoveRandAll=[];
nSpikes=zeros(1,nPatchSweeps);
nSpikesRand=zeros(1,nPatchSweeps);
allSpikeTimes=cell(1,nPatchSweeps);
allSpikeTimesS=cell(1,nPatchSweeps);
for n=1:nPatchSweeps
allSpikeTimes{n}=find(Msort(n,1,:)==1); %in sampling frequency units
allSpikeTimesS{n}=find(MsortS(n,1,:)==1); %in sampling frequency units
if ~isempty(allSpikeTimes{n})
tmpSpikes=allSpikeTimes{n};
nSpikes(n)=numel(tmpSpikes);
placesInMp=bsxfun(@plus,tmpSpikes,(1-preSpikeSamples):postSpikeSamples);
spikesToRemove=sum(placesInMp<=0,2) | sum(placesInMp>nPatchSamples,2);
placesInMp(spikesToRemove,:)=[]; %delete spikes that would have results in cut patch traces
nSpikes(n)=size(placesInMp,1);
Mtmp=reshape(Mp(n,placesInMp),nSpikes(n),winSpikeSamples);
Mall=[Mall;Mtmp];
spikesToRemoveAll=[spikesToRemoveAll;spikesToRemove];
%sweepSTA(j,:)=mean(Mtmp);
%old shuffling
%{
tmpIntervals=[ceil(allSpikeTimes{n}(1)/2);diff(allSpikeTimes{n})]; %the first spike time is divided by half since otherwise the last spike will alway be in the same place
randPermSequence=randperm(numel(tmpIntervals));
allSpikeTimesS{n}= cumsum(tmpIntervals(randPermSequence));
MShuffled(n,allSpikeTimesS{n})=1;
%}
tmpSpikes=allSpikeTimesS{n};
nSpikesRand(n)=numel(tmpSpikes);
placesInMp=bsxfun(@plus,tmpSpikes,(1-preSpikeSamples):postSpikeSamples);
spikesToRemoveRand=sum(placesInMp<=0,2) | sum(placesInMp>nPatchSamples,2);
placesInMp(spikesToRemoveRand,:)=[];
nSpikesRand(n)=size(placesInMp,1);
Mtmp=reshape(Mp(n,placesInMp),nSpikesRand(n),winSpikeSamples);
MallRand=[MallRand;Mtmp];
spikesToRemoveRandAll=[spikesToRemoveRandAll;spikesToRemoveRand];
end
end %for over sweeps
%Mall=Mall-Mall(:,preSpikeSamples)*ones(1,preSpikeSamples+postSpikeSamples); %set zero at the spike crossing point
%MallRand=MallRand-MallRand(:,preSpikeSamples)*ones(1,preSpikeSamples+postSpikeSamples); %set zero at the spike crossing point
%calculate means
if ~isempty(Mall)
mMall{k}(m,:)=mean(Mall,1);
mMallRand{k}(m,:)=mean(MallRand,1);
nSpikesAll{k}(m,:)=nSpikes;
nSpikesRandAll{k}(m,:)=nSpikesRand;
end
if calculateFullTraces
triggeredData.M{i,j,k,m}=Mall;
triggeredData.MR{i,j,k,m}=MallRand;
triggeredData.spikesToRemoveAll{i,j,k,m}=spikesToRemoveAll;
triggeredData.spikesToRemoveRandAll{i,j,k,m}=spikesToRemoveRandAll;
end
%plotShifted(((1-preSpikeSamples):postSpikeSamples)/20,Mall','verticalShift',30);line([0 0],ylim,'color',[0.8 0.8 0.8]);
%plot(((1-preSpikeSamples):postSpikeSamples)/20,mean(Mall));hold on;plot(((1-preSpikeSamples):postSpikeSamples)/20,mean(MallRand),'r');
end %for loop over Neurons (m=neuorns in series)
end%for loop over series (k)
%save to files only for large scale analysis
if nargin<=3
save(obj.Par.(Exp{i}).(Neurons{j}).seriesSortingName,'t','ic','Triggers','-v7.3');
save(obj.Par.(Exp{i}).(Neurons{j}).MEATriggeredPatchFilename,'mMall','mMallRand','nSpikesAll','nSpikesRandAll','-v7.3','-mat');
end
end %if condition to overwrite data
end %for loop over Neurons (j)
end %for loop over recordings (i)
catch errorMsg
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['An error occured while running runMEATriggeredAvgOnPatchData on ' getenv('COMPUTERNAME') ' session ' num2str(i) '/' num2str(nExp)],errorMsg.getReport);
end
rethrow(errorMsg);
end
if nargin==1 %only for long-term analysis
if obj.sendEMailMessages
sendMailViaGmail(obj.email4Messages,['runMEATriggeredAvgOnPatchData completed on ' getenv('COMPUTERNAME') ,' session success: ' num2str(nExp)]);
end
end
end %MEATriggeredAvgOnPatchData
function [obj,MEASTA]=checkTriggers(obj,Exp,Neurons,Series)
displayOnlySeriesWithTriggerIssues=true;
if nargin<2
Exp=obj.Par.experimentNames;
end
nExp=numel(Exp);
for i=1:nExp %go over all experiments
disp(['******Checking triggers for exp: ' Exp{i}]);
if nargin<3
Neurons=fields(obj.Par.(Exp{i}));
end
for j=1:numel(Neurons) %go over all Neurons within an experiment
disp(['Analyzing ' Neurons{j}]);
%check if to only analyze specific series, else analyzes all series
MEAfileNames=obj.Par.(Exp{i}).(Neurons{j}).fullMEAFilenames;
patchSeries=obj.Par.(Exp{i}).(Neurons{j}).Series;
nSeries=numel(patchSeries);
if nargin<4
pSeries=1:nSeries;
else
[~,pSeries]=intersect(cell2mat(obj.Par.(Exp{i}).(Neurons{j}).Series),Series);
nSeries=numel(pSeries);
end
for k=1:nSeries %go over series in patch data
ind=pSeries(k);
disp(['Analyzing Series ' num2str(patchSeries{ind})]);
patchUpsampledFilename=obj.Par.(Exp{i}).(Neurons{j}).patchUpsampledFilename{ind};
load(patchUpsampledFilename,'nPatchSweeps');
nTriggers=zeros(1,numel(MEAfileNames{ind}));
for m=1:numel(MEAfileNames{ind})
if ispc
MCR=MCRackRecording(MEAfileNames{ind}{m});
elseif isunix
MCR=MCRackRecordingNeuroshare(MEAfileNames{ind}{m});
end
MCR.includeDigitalDataInTriggers=true;
Triggers=MCR.getTrigger;
if isempty(Triggers)
nTriggers(m)=0;
else
nTriggers(m)=numel(Triggers{1});
end
end
if nPatchSweeps==sum(nTriggers)
if displayOnlySeriesWithTriggerIssues
disp('Series ok!');
end
else
disp('Incompatible trigger numbers!!!');
disp(['Series has ' num2str(nPatchSweeps) ' sweeps']);
disp('Mcd files are:');
disp(MEAfileNames{ind}');
disp('Corresponding trigger numbers are:');
disp(nTriggers');
end
end
end
end
end
function obj=saveParametersToFile(obj,fileName)
if nargin==1
name='parameterFile';
pathstr=cd;
else
[pathstr,name,ext] = fileparts(fileName);
end
Par=obj.Par;
if isempty(pathstr)
pathstr=cd;
end
if ispc
saveFilename=[pathstr filesep name '.mat'];
elseif isunix
saveFilename=[pathstr filesep name '_linux.mat'];
end
save(saveFilename,'Par','-v7.3','-mat');
obj.ParFilename=saveFilename;
end
function obj=loadParametersFromFile(obj,fileName)
if nargin==1
name='parameterFile';
pathstr=cd;
else
[pathstr,name,ext] = fileparts(fileName);
end
if isempty(pathstr)
pathstr=cd;
end
if ispc
loadFilename=[pathstr filesep name '.mat'];
elseif isunix
loadFilename=[pathstr filesep name '_linux.mat'];
end
load(loadFilename);
obj.Par=Par;
obj.ParFilename=loadFilename;
end
end %Methods
end %Object