Permalink
Cannot retrieve contributors at this time
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.asv
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
1575 lines (1316 sloc)
87.2 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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'};%,'ID',};%'Experiment'}; | |
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=20;%ms | |
Par.all.postMEASpikeOnPatch=100;%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 collect 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 | |
% if ~isnan([Par.(Exp{i}).(Neurons{j}).MEAfile{1}]) % Added by Lorenz, only if there's a MEA file specified, do the following lines | |
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 | |
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); | |
pSeriesNoNaN = find(cell2mat(obj.Par.(exp).(neuron).Series(~isnan(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)); | |
% channelNumbers=sort([34 35 36 32 31 28 27 26 30 1 24 23 22 18 2 20 19 14 10 6 17 16 13 9 5]); | |
% | |
% channelNumbers=sort([44 52 56 66 70 40 48 60 62 78 36 32 31 61 90 26 30 1 91 92 22 18 2 120 108]); | |
channelNumbers=1:120; | |
% channelNumbers=sort([33 34 35 36 32 31 61 90 86 87 88 89]); % one line of electrodes % | |
numberOfSpikes=1:280; | |
% a=squeeze(mean(MAll(channelNumbers,numberOfSpikes,:),2)); | |
% for i = 1 : length(channelNumbers);b(i,:)=a(i,:)-median(a(i,:));end | |
% [hPlot, scale]=activityTracePhysicalSpacePlot(h,channelNumbers',b,En,'DrawElectrodeNumbers',0); | |
% % % xlim([1.5 6.5]);ylim([2.5 7.5]) | |
% xlim([3.5 8.5]);ylim([4.5 9.5]) | |
% | |
allSpikes=cumsum(nSpikes); | |
spikeIndex=[[1;allSpikes(1:end-1)'+1] allSpikes(1:end)']; | |
clear allSpikes | |
for i = 1 :length(nSpikes) | |
allSpikes{i}=spikeIndex(i,1):spikeIndex(i,2); | |
end | |
seriesIndex= [1 3:5 7:8 10 12 13 15 17 19 21:22 24 26 28 29 30 32:40 42:50]; %[1:39 42 46:50] | |
allSpikes=allSpikes(seriesIndex); | |
allSpikes=cell2mat(allSpikes); | |
MAll1=MAll(:,allSpikes,:); | |
for i = 1 : size(MAll,2) | |
if max(MAll(1,i,:))>50 || min(MAll(1,i,:))<-50 | |
epileptic(i) = true; | |
end | |
end | |
MAll1=MAll; | |
MAll1(:,find(epileptic),:)=[]; | |
MAll1=bsxfun(@minus,MAll1,mean(mean(MAll1,2),3)); | |
[hPlot, scale]=activityTracePhysicalSpacePlot(h,channelNumbers',squeeze(mean(MAll1(1:end-1,:,:),2)),En,'DrawElectrodeNumbers',0); | |
% [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=plotPatchTriggeredMEASingleElectrode(obj,h,exp,neuron,series,preSpike_ms,postSpike_ms,electrodeNr,ylimits) | |
ylimits=[-30 5]; | |
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); | |
pSeriesNoNaN = find(cell2mat(obj.Par.(exp).(neuron).Series(~isnan(cell2mat(obj.Par.(exp).(neuron).Series))))==series); | |
set(get(h,'Parent'),'position',[466 709 668 321]); | |
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 | |
% numberOfSpikes=1:280; | |
% a=squeeze(mean(MAll(channelNumbers,numberOfSpikes,:),2)); | |
% for i = 1 : 120;b(i,:)=a(i,:)-median(a(i,:));end | |
% | |
% [hPlot, scale]=activityTracePhysicalSpacePlot(h,channelNumbers',b,En,'DrawElectrodeNumbers',0); | |
% % xlim([1.5 6.5]);ylim([2.5 7.5]) | |
% xlim([3.5 8.5]);ylim([4.5 9.5]) | |
onlyFirstSpikes=[1 cumsum(nSpikes(1:end-1))+1]; | |
secondSpikes=onlyFirstSpikes+1; | |
onlyLastSpikes=[cumsum(nSpikes(1:end))]; | |
time=-preSpike_ms+(1000/Fs):1000/Fs:postSpike_ms-40; | |
% plot(h,time,squeeze(b(electrodeNr,1:length(time))),'-k') | |
plot(h,time,squeeze(mean(MAll(electrodeNr,secondSpikes,1:(length(time))),2))','-k') | |
% plot(h,time,squeeze(mean(MAll(electrodeNr,:,1:(length(time))),2))','-k') | |
xlabel('ms');ylabel('uV'); | |
ylim(ylimits) | |
title(h,['Voltage trace, electrodeNr ' num2str(electrodeNr) ', nr of spikes' num2str(size(MAll,2))]) | |
end %plotPatchTriggeredMEASingleElectrode | |
function obj=plotPatchTriggeredMEAAndPatchSingleElectrode(obj,h,exp,neuron,series,preSpike_ms,postSpike_ms,electrodeNr,plotSpikeNr,ylimitsMea,ylimitsPatch) | |
if nargin < 9 | |
plotSpikeNr=1; | |
end | |
if nargin < 11 | |
ylimitsMEA=[-30 5]; | |
ylimitsPatch=[-60 50]; | |
end | |
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); | |
pSeriesNoNaN = find(cell2mat(obj.Par.(exp).(neuron).Series(~isnan(cell2mat(obj.Par.(exp).(neuron).Series))))==series); | |
% load patch data | |
load(cell2mat(obj.Par.(exp).(neuron).patchUpsampledFilename(pSeries))); | |
load(cell2mat(obj.Par.(exp).(neuron).patchTriggeredMEAFilename(pSeries))); | |
% set figure and add second panel | |
set(get(h,'Parent'),'position',[466 325 670 705]); | |
% hSecond=figure('position',[466 325 670 705]); | |
% hLast = figure('position',[466 325 670 705]); | |
h1=subplot(2,1,1); | |
h2=subplot(2,1,2); | |
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 | |
% numberOfSpikes=1:280; | |
% a=squeeze(mean(MAll(channelNumbers,numberOfSpikes,:),2)); | |
% for i = 1 : 120;b(i,:)=a(i,:)-median(a(i,:));end | |
% | |
% [hPlot, scale]=activityTracePhysicalSpacePlot(h,channelNumbers',b,En,'DrawElectrodeNumbers',0); | |
% % xlim([1.5 6.5]);ylim([2.5 7.5]) | |
% xlim([3.5 8.5]);ylim([4.5 9.5]) | |
if strmatch(plotSpikeNr,'last') | |
specificSpikes=[cumsum(nSpikes(1:end))]; | |
else | |
specificSpikes=[1 cumsum(nSpikes(1:end-1))+1]; | |
specificSpikes=specificSpikes + plotSpikeNr-1; | |
end | |
time=-preSpike_ms+(1000/Fs):1000/Fs:postSpike_ms; | |
% plot(h,time,squeeze(b(electrodeNr,1:length(time))),'-k') | |
counter=0; | |
for k = 1 : length(pSpikes) | |
counter=counter+1; | |
if strmatch(plotSpikeNr,'last') | |
nrOfSpikesNeeded=1; | |
else | |
nrOfSpikesNeeded=plotSpikeNr; | |
end | |
if length(pSpikes{k})>nrOfSpikesNeeded | |
if strmatch(plotSpikeNr,'last') | |
pSpikes{k}((pSpikes{k}/Fs)>4)=[]; | |
allSpecificSpikeShapes(counter,:) = Mp(k,pSpikes{k}(end)-preSpike+1:pSpikes{k}(end)+postSpike); | |
else | |
pSpikes{k}((pSpikes{k}/Fs)<2)=[]; | |
allSpecificSpikeShapes(counter,:) = Mp(k,pSpikes{k}(plotSpikeNr)-preSpike+1:pSpikes{k}(plotSpikeNr)+postSpike); | |
end | |
else | |
counter=counter-1; | |
end | |
end | |
plot(h1,time,squeeze(mean(MAll(electrodeNr,specificSpikes,1:(length(time))),2))','-k') | |
plot(h2,time,allSpecificSpikeShapes,'-k') | |
% plot(h,time,squeeze(mean(MAll(electrodeNr,:,1:(length(time))),2))','-k') | |
xlabel('ms');ylabel('uV'); | |
ylim(h1,ylimitsMEA) | |
ylim(h2,ylimitsPatch) | |
title(h1,['Voltage trace, electrodeNr ' num2str(electrodeNr) ', nr of spikes' num2str(size(MAll,2)) ' only spikes nr ' num2str(plotSpikeNr)]) | |
end %plotPatchTriggeredMEASingleElectrode | |
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); | |
pSeriesNoNaN = find(cell2mat(obj.Par.(exp).(neuron).Series(~isnan(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{pSeriesNoNaN},2); | |
else | |
disp(['Cant plot - non existing file: ' obj.Par.(exp).(neuron).patchTriggeredMEAFilename(series) 'series ' num2str(series)]); | |
return; | |
end | |
load(['layout_' obj.Par.(exp).(neuron).mea{pSeries}]); | |
%removes these lines after rerunnnig analysis | |
pCh=(nSpikes>=minNumberOfSpikes2Display); | |
channelNumbers=ic{pSeriesNoNaN}(1,pCh)'; | |
neuNumbers=ic{pSeriesNoNaN}(2,pCh); | |
STA=mMall{pSeriesNoNaN}(pCh,pInData); | |
STARand=mMallRand{pSeriesNoNaN}(pCh,pInData); | |
nSpikes=nSpikes(pCh); | |
if isempty(STA) | |
disp(['All neurons do not have enough spikes in: ' num2str(exp) '-' num2str(neuron) '_series ' num2str(series) ' ,lower the min threshold to get results']); | |
return; | |
end | |
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(pSeries)==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); | |
GS=gridSorter(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); | |
GS.sortingDir = sortingDir; | |
GS.clusteringMinSpikesTotal = 2; | |
GS.clusteringMinNSpikesCluster = 1; | |
GS.localGridSize = gridSize; | |
GS.localGridExt = gridExt; | |
GS.fastPrinting = fastPrinting; | |
GS.overwriteSpikeExtraction = overwriteSpikeExtraction; | |
GS.overwriteFeatureExtraction = overwriteFeatureExtraction; | |
GS.overwriteClustering = overwriteClustering; | |
GS.overwriteMerging = overwriteMerging; | |
GS.overwriteFitting = overwriteFitting; | |
GS.overwritePostProcessignAnalysis = overwritePostProcessignAnalysis; | |
if isempty(channels) | |
% Nothing else to set | |
% 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); | |
GS.selectedChannelSubset=channelNumbers; | |
end | |
% Start sorting | |
GS=GS.runSorting; | |
Triggers=MCR.getTrigger; | |
load(GS.sortingFileNames.fittingFile); | |
load(GS.sortingFileNames.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} '. Extracting patch data.']); | |
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{ind}) | |
load(tmpPatchDataName{ind},'data'); %it is assumned that the same Neuron will always belong to the same patch recording | |
loadedPatchFile=tmpPatchDataName{ind}; | |
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=5; | |
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} ' Extracting patch-triggered MEA data.']); | |
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{ind}) | |
load(tmpPatchDataName{ind},'data'); %it is assumned that the same Neuron will always belong to the same patch recording | |
loadedPatchFile=tmpPatchDataName{ind}; | |
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); | |
spikeLocalSmooth=csaps(1:spikePeakDetectionWinSamples,spikeLocal',0.02,1:spikePeakDetectionWinSamples,ones(1,spikePeakDetectionWinSamples))'; | |
[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 | |
%figure('Position', [50 50 1200 800]);plot(Mp(m,:));hold on;plot(patchSpikesLogical*3,'r');plot(med+obj.Par.all.medFiltDev*medDev,'g');plot(pSpikes{m},Mp(m,pSpikes{m}),'xk');xlim([0.4 0.8]*10^5) | |
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; | |
MEASTA.pSpikes=pSpikes; | |
MEASTA.nSpikes=nSpikes; | |
end | |
elseif strcmp(clamp{ind},'VC') | |
disp('Series is in voltage clamp mode. No extraction of spikes possible, therefore nothing to extract. No calculations performed!'); | |
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} '. Extracting MEA-spike-triggered data from patch recording.']); | |
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 | |
if nSeries>0 | |
%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 | |
%} | |
% LORENZ CHANGED LINE: | |
% Triggers{k}=TriggersAll{1}(startTrigSweep(ind):endTrigSweep(ind))-delay2Trig{ind}*obj.resamplingFactor;%get triggers and | |
Triggers{k}=TriggersAll{1}(startTrigSweep(k):endTrigSweep(k))-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 nSeries>0 %number of valid series > 0 | |
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 |