Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % NOTES
- % SpikeMat (ts x channels )
- % Creates and populates UnitQualityMat with Loc,Name,Ch,Unit Info
- % v2: add waveforms to the unitInfoMat
- % nhpPath = goto_nhp();
- % dataPath = fullfile(nhpPath, 'data');
- % goto_nhp; cd('data');
- %% Edit ------------------------------------------------------------------------
- cd('C:\Users\rocke\Desktop\nhp-project')
- dbstop if error
- if nargin == 0
- fileDirs = {
- './data/2016-05-18_jones'
- };
- end
- fs = 1000; %Hz, resolution of spike proccess
- overwriteBool = false;
- allAlignments = {'recordingStart', 'condTone', 'pureTone', 'allTone', 'puff'};
- % allAlignments = {'recordingStart', 'DBSonset', 'DBSoffset', 'condTone', 'pureTone', 'allTone', 'puff'}
- minQuality = .8;
- % monkeyname = 'mr. jones';
- % secsBeforeDataStart = 50;
- % secsPostDataStart = 130;
- debugBool = 0;
- %% -----------------------------------------------------------------------------
- %% DEV TODOs
- % ...
- %% parse options
- if nargin
- % ensure input is cell array
- if ~iscell(fileDirs)
- fileDirs = {fileDirs};
- end
- end
- % ensure allAlignments is cell array
- if ~iscell(allAlignments)
- allAlignments = {allAlignments};
- end
- nAlignments = length(allAlignments);
- nDirs = length(fileDirs);
- for iDir = 1:nDirs
- fileDir = fileDirs{iDir};
- fprintf('Dir (%i/%i): %s \n', iDir, nDirs, fileDir);
- % check for existing files
- saveFilesExist = nan(nAlignments,1);
- for iAlignment = 1:nAlignments
- thisAlignment = allAlignments{iAlignment};
- saveFilename = sprintf('spikeDataCombined_%s.mat', thisAlignment); % filename based on alignment
- saveFilepath = fullfile(fileDir, saveFilename);
- saveFilesExist(iAlignment) = exist(saveFilepath, 'file');
- clear thisAlignment saveFilename saveFilepath
- end
- if all(saveFilesExist) && ~overwriteBool
- fprintf(' All save files exist, skipping dir... \n');
- % continue
- end
- % Get file names
- files = lscell(fileDir);
- dataFiles = files(contains(files, '.plx'));
- nFiles = length(dataFiles);
- % skip if no files
- if isempty(dataFiles)
- % continue
- end
- %% unitInfoMatCombined
- unitInfoMatCombined = load(fullfile(fileDir, 'unitInfoMatCombined.mat'));
- unitInfoMatCombined = unitInfoMatCombined.UnitInfoMatCombined;
- %% check number of units
- for iFile = 1:nFiles
- file = dataFiles{iFile};
- filePathIn = fullfile(fileDir, file);
- spikeCounts{iFile} = plx_info(filePathIn, 1);
- unitsPerCh{iFile} = sum(logical(spikeCounts{iFile}));
- totalUnitsPerFile(iFile) = sum(unitsPerCh{iFile});
- end
- if sum(totalUnitsPerFile) ~= length(unitInfoMatCombined)
- warning('Number of units found does not match those known from sessionInfo file. Skipping dir...')
- % continue
- end
- %% filter units
- qualityUnits = ([unitInfoMatCombined.Quality] >= minQuality);
- qualityUnits = qualityUnits & [unitInfoMatCombined.Isolated]; % check isolated==true
- % nQualityUnits = sum(qualityUnits);
- assert(nFiles == 2)
- % split into 2 files
- qualityUnitsInFileLogical{1} = qualityUnits(1:totalUnitsPerFile(1));
- qualityUnitsInFileLogical{2} = qualityUnits(totalUnitsPerFile(1)+1:end);
- % update unitInfoMatCombined
- unitInfoMatCombined(~qualityUnits) = [];
- %% load plexon spike data (*.plx)
- fprintf(' Loading plexon spike data... \n')
- % Open files
- spikeTimesSec = cell(nFiles,1);
- spikeCounts = cell(nFiles,1);
- for iFile = 1:nFiles
- % File Path
- file = dataFiles{iFile};
- filePathIn = fullfile(fileDir, file);
- [~,fileName] = fileparts(filePathIn);
- % File Print
- fprintf(' File (%i/%i): %s \n', iFile, nFiles, fileName);
- qualityUnitsInFile{iFile} = find(qualityUnitsInFileLogical{iFile});
- % load data
- verboseBool = false;
- [spikeTimesSec{iFile}, ~, ~, durationSec, plxFs] = nhpPlxDataExtract(filePathIn, verboseBool, debugBool, qualityUnitsInFile{iFile});
- % spikeTimesSec{file}{unit}
- % Note: assumes last ch loaded by plx is empty/extra
- end
- fprintf(' Done loading plexon spike data. \n');
- %% load sessionInfo
- sessionInfoFile = files{contains(files, 'sessionInfo')};
- sessionInfoPath = fullfile(fileDir, sessionInfoFile);
- sessionInfoData = load(sessionInfoPath);
- sessionInfo = sessionInfoData.sessionInfo;
- %% create binary spike process
- fprintf(' Creating binary spike process... \n');
- for thisAlignment = allAlignments(:)'
- thisAlignment = thisAlignment{1};
- fprintf(' Alignment: %s \n', thisAlignment);
- % start metadata
- metadata.alignment = thisAlignment;
- metadata.unitInfoMatCombined = unitInfoMatCombined;
- metadata.sessionInfo = sessionInfo;
- % check for existing file
- saveFilename = sprintf('spikeDataCombined_%s.mat', thisAlignment); % filename based on alignment
- saveFilepath = fullfile(fileDir, saveFilename);
- if exist(saveFilepath, 'file') && ~overwriteBool
- fprintf(' Save file exists, skipping alignment %s (%s)\n', thisAlignment, saveFilepath);
- % continue
- end
- spkDataCombined = cell(nFiles, 1);
- for iFile = 1:nFiles
- thisFileSpikeTimesSec = spikeTimesSec{iFile};
- % nonEmptyChs = ~cellfun(@isempty, thisFileSpikeTimesSec);
- % unitsPerCh = sum(logical(spikeCounts{iFile}));
- % nTotalUnits = sum(unitsPerCh{iFile});
- nTotalQualityUnits = length(qualityUnitsInFile{iFile});
- % ignore the empty last extra channel erroneously added by readall plx
- % nonEmptyChs(end) = 1;
- switch thisAlignment
- case 'recordingStart'
- unitDim = 2;
- nTimeInd = ceil(durationSec * fs);
- spkDataCombined{iFile} = zeros(nTimeInd, nTotalQualityUnits, 'single'); % time by units
- for iUnit = 1:nTotalQualityUnits
- thisUnitSpikeTimesSec = thisFileSpikeTimesSec{iUnit};
- % convert times to indicies
- thisUnitSpikeInds = round(thisUnitSpikeTimesSec * fs);
- spkDataCombined{iFile}(thisUnitSpikeInds, iUnit) = 1;
- end
- clear iUnit
- case {'condTone', 'pureTone', 'allTone', 'puff'}
- % cond tone start min about 1 sec away from puff
- % pure tone start min about 3 sec away from puff
- % 2 tones start min about 4 sec away from each other
- % tone: 500 ms
- % puff: ~106 ms (105-110 ms skewed right)
- % convert trial start time stamps to unts of ms
- puffTimesSec = sessionInfo.CondPuffTone.PuffTimes / plxFs;
- condToneTimesSec = sessionInfo.CondPuffTone.condSoundOnsetTimesPts / plxFs;
- pureToneTimeSec = sessionInfo.CondPuffTone.pureToneOnsetTimesPts / plxFs;
- unitDim = 3;
- switch thisAlignment
- case 'condTone'
- % 500 ms before start, 2000 ms after start
- % gets 500 ms before tone, tone, puff, and ~1400 ms after puff
- preEvntTimeLenInMs = 500; %ms
- postEvntTimeLenInMs = 2000; %ms
- trialEvntsSec = condToneTimesSec;
- spkDataCombined{iFile} = makeTrialData(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- metadata.trials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- case 'pureTone'
- % 500 ms before start, 2000 ms after start
- % gets 500 ms before tone, tone, and 1500 ms after tone, but no puff
- preEvntTimeLenInMs = 500;
- postEvntTimeLenInMs = 2000;
- trialEvntsSec = pureToneTimeSec;
- spkDataCombined{iFile} = makeTrialData(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- metadata.trials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- case 'allTone'
- % 500 ms before start, 900 ms after start
- % gets 500 ms before tone, tone, and 400 ms after tone, but no puff
- preEvntTimeLenInMs = 500; %ms
- postEvntTimeLenInMs = 900; %ms
- trialEvntsSec = sort([condToneTimesSec(:); pureToneTimeSec(:)]);
- spkDataCombined{iFile} = makeTrialData(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- metadata.trials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- metadata.pureTrials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, pureToneTimeSec);
- metadata.condTrials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, condToneTimesSec);
- case 'puff'
- % 500 ms before start, 1100 ms after start
- % gets 500 ms before puff, puff and ~1000 ms puff tone, but no tones
- preEvntTimeLenInMs = 500; %ms
- postEvntTimeLenInMs = 1100; %ms
- trialEvntsSec = puffTimesSec;
- spkDataCombined{iFile} = makeTrialData(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- metadata.trials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec);
- end % 2nd switch alignment
- case {'DBSonset', 'DBSoffset'}
- unitDim = 3;
- % FIXME
- % nDBS = length(sessionInfo.DBS);
- % nBins = (secsBeforeDataStart+secsPostDataStart)*1000;
- % spkData = zeros(nBins,totalUnits,nDBS);
- %
- % % make spkdata equivalent size with interneuron/unit quality/ unsorted
- % tempTrialSpks = []; %binSpikeCountsTrial = zeros(nBins,1);
- % unitNumCrossRef = [];
- % chNameUnitCrossRef = {};
- % unitWaveforms = {};
- %
- % % MAYBE USE HIST INSTEAD???
- % for d = 1:nDBS
- % if contains(alignment,'DBSonset','IgnoreCase',true)
- % trialZeroMsec = round(sessionInfo.DBS(d).onsetOffsetTimePts(1,1)/30);
- % % %PuffsAdjustedTimes = [LOCtimePts - dataZeroPt + secsBeforeDataStart*fs]; % time of LOC (in drugStart dataMatrix) from the onset of the data clip
- % % put in cell array
- % elseif contains(alignment,'DBSoffset','IgnoreCase',true)
- % trialZeroMsec = round(sessionInfo.DBS(d).onsetOffsetTimePts(1,2)/30);
- % end
- % %Trial Start in Msec Pts
- % dataStartPts = trialZeroMsec - secsBeforeDataStart*fs; % points (msec for NS2) before (drug) start to take data
- % dataEndPts = trialZeroMsec + secsPostDataStart*fs;
- % unitCount = 0;
- % for ch = 1:size(allts,2)-1 % %subtract the empty extra channel erroneously added by readall plx
- % for u=1:maxUnitsPerCh(ch)
- % unitCount = unitCount+1;
- % unitNumCrossRef(unitCount,:) = [ch u];
- % chNameUnitCrossRef{unitCount} = channel_names{ch};
- % binSpikeCountsTrial = zeros(nBins,1);
- % tempTrialSpks = allts{u,ch}(allts{u,ch}> (dataStartPts/fs) & allts{u,ch}< (dataEndPts/fs));
- % if isempty(tempTrialSpks)
- % else
- % tempAdjSpkTimes = round(tempTrialSpks*fs)-dataStartPts;
- % if tempAdjSpkTimes(1) == 0 % if first spike happens at time 0 of trial (toss it)
- % tempAdjSpkTimes = tempAdjSpkTimes(2:end); %toss first spike, cant index @ 0
- % end
- % unitWaveforms{unitCount} = allwf{u,ch};
- % binSpikeCountsTrial([tempAdjSpkTimes])=1; %spikes adjusted to bins
- % spkData(:,unitCount,d) = binSpikeCountsTrial;
- % % see if short ISIs, warn how many
- % shortISIind = find(diff(tempTrialSpks)<.001);
- % if any(shortISIind)
- % warning('DBS run %d: %d short ISI spikes for unit %d ',d,length(shortISIind),unitCount)
- % % for total n spikes do nnShortISIS =1 for each count
- % end
- % end
- % end
- % end
- % end
- % end % case {'DBSonset', 'DBSoffset'}
- end % switch alignment
- %% FIXME
- % %% indexing into the unit data matrix
- % chNamesSTR = string(channel_names);
- % elec1ChInd = contains(chNamesSTR,'elec1');
- % elec2ChInd = contains(chNamesSTR,'elec2');
- %
- % if strcmpi(arrayLocation,'pfc')%elec 1- :FEF, 2- vlPFC
- % sprintf('making spikeMat for pfc: %s',monkeyname)
- % FEFchInd = elec1ChInd;
- % vlPFCchInd = elec2ChInd;
- % FEFchIndNum = find(FEFchInd ==1);
- % FEFUnits = ismember(unitNumCrossRef(:,1),FEFchIndNum);
- % vlPFCchIndNum = find(vlPFCchInd ==1);
- % vlPFCUnits = ismember(unitNumCrossRef(:,1),vlPFCchIndNum);
- % elseif strcmpi(arrayLocation,'ppcstg')
- % sprintf('making spikeMat for ppcstg: %s',monkeyname)
- % %mr. jones: 1- stg; 2-ppc; mary 2- stg, 1- ppc
- % STGchInd = elec1ChInd;
- % PPCchInd = elec2ChInd;
- % if contains(monkeyname,'mary','IgnoreCase',true) %flip the elec assignments
- % STGchInd = elec2ChInd;
- % PPCchInd = elec1ChInd;
- % end
- % PPCchIndNum = find(PPCchInd ==1);
- % PPCUnits = ismember(unitNumCrossRef(:,1),PPCchIndNum);
- % STGchIndNum = find(STGchInd ==1);
- % STGUnits = ismember(unitNumCrossRef(:,1),STGchIndNum);
- % end
- %
- % %% Make unitInfoMat (each entry is unit in this spike matrix)
- % % populates: Loc, ElecName,Ch,UnitNum (in its ns6/.dat file),... TimeRangeSecs (0 = unstable/in and out sort), Waveforms
- % unitInfoMat(1:totalUnits) = struct('Location',[],'ElectrodeName',[],'Channel',[],'UnitNum',[],'Quality',[],'Interneuron',[],'Isolated',[],'TimeRangeSecs',[],'Waveforms',[]);
- %
- % if contains(arrayLocation,'pfc','IgnoreCase',true)
- % [L1{1:length(FEFUnits)}] = deal('FEF');
- % [unitInfoMat(FEFUnits).Location] = L1{:};
- % [L2{1:length(vlPFCUnits)}] = deal('vlPFC');
- % [unitInfoMat(vlPFCUnits).Location] = L2{:};
- % else
- % [L1{1:length(PPCUnits)}] = deal('PPC');
- % [unitInfoMat(PPCUnits).Location] = L1{:};
- % [L2{1:length(STGUnits)}] = deal('STG');
- % [unitInfoMat(STGUnits).Location] = L2{:};
- % end
- % [unitInfoMat(:).ElectrodeName] = chNameUnitCrossRef{:}; %true for both array sets
- % unitNumXrefCell = num2cell(unitNumCrossRef);
- % [unitInfoMat(:).Channel] = unitNumXrefCell{:,1};% clarify unit num 1 = unsorted
- % [unitInfoMat(:).UnitNum] = unitNumXrefCell{:,2};
- % [unitInfoMat(:).Waveforms] = unitWaveforms{:};
- end % files
- % concat units from both files into numeric array from cell array
- fileDataLengths = cellfun(@(x) size(x,1), spkDataCombined);
- % trim data to min length if not equal lengths
- if ~all(~diff(fileDataLengths))
- minLen = min(fileDataLengths);
- spkDataCombined = cellfun(@(x) x(1:minLen, :, :), spkDataCombined, 'Uni', false);
- end
- % combine files
- spkDataCombined = cat(unitDim, spkDataCombined{:}); %#ok<NASGU>
- % Save Spike Data (not binned, well only 1msec bin)
- % save(saveFilepath, 'spkDataCombined','metadata','-v7.3');
- fprintf(' Saved: %s \n', saveFilepath);
- clear spkDataCombined metadata
- end % allAlignments
- end % dirs
- fprintf('\nDone! \n\n')
- %% Nested fn
- function spkDataCombined = makeTrialData(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec)
- % parse inputs
- preEvntTimeIndLen = round(preEvntTimeLenInMs * fs/1000);
- postEvntTimeIndLen = round(postEvntTimeLenInMs * fs/1000);
- trialEvntsInd = round(trialEvntsSec * fs);
- nTrials = length(trialEvntsInd);
- % preallocate
- nTimeInd = preEvntTimeIndLen + postEvntTimeIndLen;
- spkDataCombined = zeros(nTimeInd, nTrials, nTotalQualityUnits, 'single'); % time by trials by units
- for iUnit = 1:nTotalQualityUnits
- thisUnitSpikeTimesSec = thisFileSpikeTimesSec{iUnit};
- % convert times to indicies
- thisUnitSpikeInds = round(thisUnitSpikeTimesSec * fs);
- for iTrial = 1:nTrials
- thisTrialEvntInd = trialEvntsInd(iTrial);
- thisTrialStartInd = (thisTrialEvntInd-preEvntTimeIndLen);
- thisTrialEndInd = (thisTrialEvntInd + postEvntTimeIndLen);
- thisTrialSpikeInds = thisUnitSpikeInds( thisUnitSpikeInds >= thisTrialStartInd & thisUnitSpikeInds <= thisTrialEndInd ) - thisTrialStartInd + 1; % +1 so doesn't start at 0
- spkDataCombined(thisTrialSpikeInds, iTrial, iUnit) = 1;
- end
- end
- clear iUnit
- end
- end % main
- %% Local fn
- function trials_Start_Evnt_End_sec = getTrialsStart_Evnt_End(preEvntTimeLenInMs, postEvntTimeLenInMs, trialEvntsSec)
- trialEvntsSec = trialEvntsSec(:); % make col vec
- % make matrix of [Start Evnt End] cols and rows of trials
- trials_Start_Evnt_End_sec = [trialEvntsSec - preEvntTimeLenInMs/1000, trialEvntsSec, trialEvntsSec + postEvntTimeLenInMs/1000];
- end
- %Finish Manual Entries
- % Quality: 0-1; Interneuron (by eye): 0/1; Isolated: 0/1;
- % TimeRangeSecs:[# #;# #]; 0 = unstable; consider? 1 = present throughout ;
- %quality of non-isolated indicates quality
- % of multiunit remaining in unsorted
- %Reorder and save spikeDataMat
- % next analyses will require:
- % c = extractfield(UnitInfoMat,'Location')
- % concatenate structs (PFC and PPCSTG) via catStruct = [unitQuality, unitQuality2];
- %find best rated neuron on single channel; must MUST take neurons with
- %time-range that ~= 0; time range of 0 means unstable
- %function to guassian kernel smooth spike train (; make kermel w/ normpdf( LOOK @ p320 in ch20 in M4Nv2
- % LATER FIGURE OUT: deal with short ISI spikes on same unit (< 1ms apart)
- %right now can incorrectly add extra spikes, need to adjust stamps before putting in binSpikeCounts
- % shortISIind = find(diff(tempTrialSpks)<.001);
- % floorTS = floor(tempTrialSpks(shortISIind)*fs)/fs;
- % ceilTS = ceil(tempTrialSpks(shortISIind+1)*fs)/fs;
- % if any(shortISIind)
- % warning('added short ISI spike to next bin for unit %d',unitCount)
- % % for total n spikes do nnShortISIS =1 for each count
- % end
- % nlostspikes = total - nnz
- % % If want actual counts per bin)
- % for b = 1:nBins
- % window = 1+(b-1): 1+b;
- % spkData(b,unitCount,n) = any(tempTrialSpks>= + binsize(;
- % spkDataMulti % dataMat w/ multiple spiikes per msec counted
- % end
- % end
- % end
- %
- % end
- %
- % binSpikeCounts = zeros(nBins,nCh,nTrials);
- % binSpikeCountsTrial = zeros(nBins,1);
- % window = [];
- % for k = 1:nTrials
- % for j = 1:nCh
- % for b = 1:nBins
- % window = 1+(b-1): 1+b;
- % window = window(window <= size(spikeMatrix,1));
- % winStarts(i) = window(1);
- % binSpikeCounts(i,j,k) = sum(spikeMatrix(window,j,k));
- % end
- % end
- % end
- % for n = 1:nDBS
- % for ch = 1:size(allts,2)-1 % %subtract the empty extra channel erroneously added by readall plx
- % for u=1:maxUnits(ch)
- % unitCount = unitCount+1;
- % spkData(:,unitCount,n) = allts(ch,u);
- % end
- % end
- % end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement