% Phase spike phase analysis (spike-time vs waveform phase by hilbert) % Execute 'A_NsxFilter_ALLchannels_v2.m' first if requires filtering of data % v2: add polarity index (Check PolarityProfiles 2) % for SW: continuous .ns6, able to choose time frame(one by one) % for SW v2: create folders & polarity index(PI) ** need validation % v2.1: polarity index=abs(polarity sum)/length ** "abs" of complex sum % => NOT the problem (changed back to without abs in this step, cuz will be caculated at polarity profile 2 eventually) % normalized to the wrong length(normalized to all spikes, not spikes in that time frame) % v2.2: check NSx&NEV sampling rate % v2_3: 2022.01.17 add 'kilo' option, uses F_openKilosort2 func./ cd before loading filtered data % v2_3b: 2022.01.18 modified a bit for filtered data from batch processing (line 71) % v2_3b2: 2022.02.07 compatible with decimated data % v2_3b3: 2022.04.07 multiple time frame & label on fig % v2_3b4: 2022.04.24 output shuffled spike phase & save a profile compilation % v2_3b5: 2022.05.16 compatible with .nev file % 2022.05.18 discard frames with too few units (below a threshold) % v2_3b6: 2022.05.22 allow filtered data with reference subtracted clear all, close all; %% Input parameters %% % LFP channel profiles LFPChannels=[10,11,129]; % can be multiple channels RefSubtracted=1; % include filtered data with Ref subtracted or not RefSubCh=[201,202,203]; % starts with 201, pseudo channel name, for indexing the Ref pair LoadFilterFile=1; % if 1, loadfilterfile for the LFP waveform kilo=1; % load kilosort data or not boss=0; % load .nev from boss %%% timeframe input at line 75 % NEV channel profiles NEVchannels=[1,1,1,1,1,1,1]; % Can be multiple channels, matrix length match NEVunits NEVunits=[6,7,8,9,10,11,12]; % Corresponding detected units of each channel, can be more than one unit in each channel % Bin spike-to-phase profiles BinNumber=36; Alltimeframes=[177,252;503,538;68,119;409,455;565,621];% in secs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% frametype=[2,2,1,1,1]; % 1: move/ 2: rest/ 3: baseline move/ 4: baseline rest HistThr=150; % if there are fewer than 150 spikes for a unit in a certain time frame, it will be saved seperately and won't be added to the total compilation %%%%%%%%% Put input parameters into structure %%%%%%%%%%%%% save temp.mat; a_InputParameters=load('temp.mat'); delete temp.mat; %% Main function %% open NEV file %% %openNEV; if kilo==1 NEV=F_openKilosort2; NEVchannels=ones(size(NEVunits)); elseif kilo==0 if boss==0 NEV=F_openNEX; elseif boss==1 NEV=openNEV; end end %% Define workspace SpikeFreqCell=cell(6,length(NEVunits)+1); s_SpikeFreqCell=cell(6,length(NEVunits)+1); A=num2cell(NEVchannels); B=num2cell(NEVunits); SpikeFreqCell(1,:)=['NEVchannels',A]; SpikeFreqCell(2,:)=['NEVchannels',B]; SpikeFreqCell(3,1)={'Phase of SU'}; SpikeFreqCell(4,1)={'SpikeHistCounts'}; SpikeFreqCell(5,1)={'PolarityCounts'}; SpikeFreqCell(6,1)={'PolarityIndex'}; s_SpikeFreqCell(1,:)=['NEVchannels',A]; s_SpikeFreqCell(2,:)=['NEVchannels',B]; s_SpikeFreqCell(3,1)={'Phase of SU'}; s_SpikeFreqCell(4,1)={'SpikeHistCounts'}; s_SpikeFreqCell(5,1)={'PolarityCounts'}; s_SpikeFreqCell(6,1)={'PolarityIndex'}; clear A B; names=["Channels";"Units";"LFPChannels";"PolarityCounts";"PolarityIndex";"Degree(-pi:pi)";"Degree(0:2pi)"]; Failed_timeframes={}; %% Load LFP file%% if LoadFilterFile==1; [A,B]=uigetfile('*.mat'); cd (B); load(A); a_InputParameters_LFPdataSource=[A,B]; clear A B; else NSx=openNSx; SourceFile=[NSx.MetaTags.Filename NSx.MetaTags.FileExt]; data=double(NSx.Data'); NSx=rmfield(NSx,'Data'); fs=NSx.MetaTags.SamplingFreq; Title1=SourceFile; end a_InputParameters.SourceFile=SourceFile; %A=find(ismember(NSx.MetaTags.ChannelID,LFPChannels))'; % modified for filtered data from batch clear i if RefSubtracted==0 A=find(ismember(filteronly,LFPChannels)); data=data(:,A); else A=find(ismember(filteronly,LFPChannels)); refA=find(ismember(str2double(withref_pseudo_ch(2:end,1)),RefSubCh))'; A=[A,refA+length(filteronly)] data=data(:,A); LFPChannels=[LFPChannels RefSubCh]; end % %% NSx SourceFile identification %% % % check if the .nev & .nsx file matched % if NSx.MetaTags.DateTime==NEV.MetaTags.DateTime; % NEV.MetaTags.Filename=NSx.MetaTags.Filename; if NEV.MetaTags.SampleRes ~= NSx.MetaTags.SamplingFreq error(sprintf('Warning: Different sampling rate, NSx:%d, NEV:%d',NSx.MetaTags.SamplingFreq,NEV.MetaTags.SampleRes)); end % if NSx.MetaTags.FileExt=='.ns6'; % else % correct the error of NEV if not processed by .ns6 file. % NEV.MetaTags.SampleRes=NSx.MetaTags.SamplingFreq; % NEV.MetaTags.DataDurationSec=NEV.MetaTags.DataDuration/NEV.MetaTags.SampleRes; % end % else % error('NEV & NSx file mismatch'); % end compilation{1,1}=SourceFile; compilation{6,1}="Shuffled"; %% Hilbert transform of LFP data %% tic HilbertLFP=hilbert(data); % hilbert transformed each column as a dataset (except only one line). shuffled_phase=HilbertLFP(randperm(length(HilbertLFP)),:); toc LegendMarker2=''; %% NEV units [sz,~]=size(Alltimeframes); for z=1:sz SpikeFreqCell=cell(6,length(NEVunits)+1); s_SpikeFreqCell=cell(6,length(NEVunits)+1); A=num2cell(NEVchannels); B=num2cell(NEVunits); SpikeFreqCell(1,:)=['NEVchannels',A]; SpikeFreqCell(2,:)=['NEVchannels',B]; SpikeFreqCell(3,1)={'Phase of SU'}; SpikeFreqCell(4,1)={'SpikeHistCounts'}; SpikeFreqCell(5,1)={'PolarityCounts'}; SpikeFreqCell(6,1)={'PolarityIndex'}; s_SpikeFreqCell(1,:)=['NEVchannels',A]; s_SpikeFreqCell(2,:)=['NEVchannels',B]; s_SpikeFreqCell(3,1)={'Phase of SU'}; s_SpikeFreqCell(4,1)={'SpikeHistCounts'}; s_SpikeFreqCell(5,1)={'PolarityCounts'}; s_SpikeFreqCell(6,1)={'PolarityIndex'}; clear A B; Failed_timeframes{z,2}=[]; timeframe=Alltimeframes(z,:); if frametype(z)==1 frametypetext='move'; elseif frametype(z)==2 frametypetext='rest'; elseif frametype(z)==3 frametypetext='baseline move'; elseif frametype(z)==4 frametypetext='baseline rest'; end dataframe=timeframe*fs; AngleEdges=-pi:2*pi/BinNumber:pi; % create Bin edges for polar system for h=1:length(NEVunits); SU_Name(h)={[LegendMarker2,' ',sprintf('Chan%d-U%d',NEVchannels(h),NEVunits(h))]}; C=find(NEV.Data.Spikes.Electrode==NEVchannels(h)); % find spikes index that match the selected ChannelID D=find(NEV.Data.Spikes.Unit(C)==NEVunits(h)); % find spikes index that match the sorted UnitID. E=NEV.Data.Spikes.TimeStamp(C(D)); % find spikes that match the sorted channel-unit. if filtertype=="decimate_fir" E=round(E/(NSx.MetaTags.SamplingFreq/dsfs)); display('**Downsampled Spikes frames'); end spikesinframe=E(E>dataframe(1)&E<dataframe(2)); %F=E/3 %%%nope, but how to transform a 30k data into 10k? % disp(size(C)); % SUpoints(h,:)=zeros(1,NEV.MetaTags.DataDuration); % plot points of SU in the timeline, step 1: fill with zeros % SUpoints(h,E)=1; % plot points of SU in the timeline %A=SpikePhase(:,:rad2deg(angle(HilbertLFP(E))); % fine [-180,180] angle each phase SpikePhase=angle(HilbertLFP(spikesinframe,:))+pi/2; % correct the -90 degree shift after hilbert transform SpikePhase=wrapToPi(SpikePhase); % wrap angle to [-pi pi] range. s_SpikePhase=angle(shuffled_phase(spikesinframe,:))+pi/2; % correct the -90 degree shift after hilbert transform s_SpikePhase=wrapToPi(s_SpikePhase); % wrap angle to [-pi pi] range. if length(spikesinframe)<HistThr Failed_timeframes{z,1}=[timeframe(1),timeframe(2)]; Failed_timeframes{z,2}=vertcat(Failed_timeframes{z,2},[NEVchannels(h),NEVunits(h)]); end SpikeFreqCell{3,h+1}=SpikePhase; %Store raw angle data into cell arrays(because the array size of each SU is different, need to use cell) for j=length(LFPChannels):-1:1 SpikeHistCounts(:,j)=histcounts(SpikePhase(:,j),AngleEdges); end SpikeFreqCell{4,h+1}=SpikeHistCounts; % save SpikeHistCount to cells PolaritySums=sum(exp(i*(SpikePhase))); % absolute vector sum of each angle (complex number) PolarityIndex=PolaritySums/length(spikesinframe); % normalized vector sum, should between [0,1] SpikeFreqCell{5,h+1}=PolaritySums; SpikeFreqCell{6,h+1}=PolarityIndex; % for shuffled data s_SpikeFreqCell{3,h+1}=s_SpikePhase; %Store raw angle data into cell arrays(because the array size of each SU is different, need to use cell) for j=length(LFPChannels):-1:1 s_SpikeHistCounts(:,j)=histcounts(s_SpikePhase(:,j),AngleEdges); end s_SpikeFreqCell{4,h+1}=s_SpikeHistCounts; % save SpikeHistCount to cells s_PolaritySums=sum(exp(i*(s_SpikePhase))); % absolute vector sum of each angle (complex number) s_PolarityIndex=s_PolaritySums/length(spikesinframe); % normalized vector sum, should between [0,1] s_SpikeFreqCell{5,h+1}=s_PolaritySums; s_SpikeFreqCell{6,h+1}=s_PolarityIndex; end [rf,~]=size(Failed_timeframes{z,2}); % Calculate Polarity count and Polarity index PolarityProfiles=zeros(5,length(NEVunits)*length(LFPChannels)); T=length(LFPChannels); for k=1:1:length(NEVunits); PolarityProfiles(1,((k-1)*T+1):k*T)=NEVchannels(k); PolarityProfiles(2,((k-1)*T+1):k*T)=NEVunits(k); PolarityProfiles(3,((k-1)*T+1):k*T)=LFPChannels; PolarityProfiles(4,((k-1)*T+1):k*T)=SpikeFreqCell{5,k+1}; % Polarity counts PolarityProfiles(5,((k-1)*T+1):k*T)=SpikeFreqCell{6,k+1}; % Polarity Index end PolarityProfile2=abs(PolarityProfiles); PolarityProfile2(6,:)=rad2deg(angle(PolarityProfiles(4,:))); % Polarity index: degree [-180,180] PolarityProfile2(7,:)=wrapTo360(PolarityProfile2(6,:)); % Polarity index: degree[0 360] PolarityProfile2=[names PolarityProfile2]; % for shuffled data s_PolarityProfiles=zeros(5,length(NEVunits)*length(LFPChannels)); T=length(LFPChannels); for k=1:1:length(NEVunits); s_PolarityProfiles(1,((k-1)*T+1):k*T)=NEVchannels(k); s_PolarityProfiles(2,((k-1)*T+1):k*T)=NEVunits(k); s_PolarityProfiles(3,((k-1)*T+1):k*T)=LFPChannels; s_PolarityProfiles(4,((k-1)*T+1):k*T)=s_SpikeFreqCell{5,k+1}; % Polarity counts s_PolarityProfiles(5,((k-1)*T+1):k*T)=s_SpikeFreqCell{6,k+1}; % Polarity Index end s_PolarityProfile2=abs(s_PolarityProfiles); s_PolarityProfile2(6,:)=rad2deg(angle(s_PolarityProfiles(4,:))); % Polarity index: degree [-180,180] s_PolarityProfile2(7,:)=wrapTo360(s_PolarityProfile2(6,:)); % Polarity index: degree[0 360] s_PolarityProfile2=[names s_PolarityProfile2]; %% Graphical illustration (polarhistogram) GetFilePath=[pwd '\'] SubPath=[sprintf('%d-%d(%s)',timeframe(1),timeframe(2),frametypetext)]; SubPath=[SubPath,'\'] if kilo==1 FullFilePath=[GetFilePath,'(kilo)SpikePhaseAnalysisTimeFrames(v2_3b6)\',SubPath]; elseif kilo==0 FullFilePath=[GetFilePath,'SpikePhaseAnalysisTimeFrames(v2_3b6)\',SubPath]; end mkdir(FullFilePath) failedfilepath=[FullFilePath,'TooFewSpikes\']; mkdir(failedfilepath) K=0; order=1; fID=1; for h=1:1:length(NEVunits); if fID<=rf if NEVchannels(h)==Failed_timeframes{z,2}(fID,1)&NEVunits(h)==Failed_timeframes{z,2}(fID,2) fID=fID+1; SpikePhase=SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(SpikePhase(:,j),BinNumber); title([sprintf('(NA)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',PolarityProfile2(5,order+1))); hold off; saveas(K,[failedfilepath,sprintf('(NA)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[failedfilepath,sprintf('(NA)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end else SpikePhase=SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(SpikePhase(:,j),BinNumber); title([sprintf('%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',PolarityProfile2(5,order+1))); hold off; saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end end else SpikePhase=SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(SpikePhase(:,j),BinNumber); title([sprintf('%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',PolarityProfile2(5,order+1))); hold off; saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[FullFilePath,sprintf('SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end end end for t=1:rf A=find(str2double(PolarityProfile2(1,:))==Failed_timeframes{z,2}(t,1)); B=find(str2double(PolarityProfile2(2,A))==Failed_timeframes{z,2}(t,2)); PolarityProfile2(:,A(B))=[]; PolarityProfiles(:,A(B)-1)=[]; A=find(cell2mat(SpikeFreqCell(1,2:end))==Failed_timeframes{z,2}(t,1))+1; B=find(cell2mat(SpikeFreqCell(2,A))==Failed_timeframes{z,2}(t,2)); SpikeFreqCell(:,A(B))=[]; end filename=sprintf('SpikePhaseAnalysis_%d-%ds_%s.mat',timeframe(1),timeframe(2),frametypetext); save([FullFilePath,filename],'SourceFile','SpikeFreqCell','PolarityProfiles','PolarityProfile2'); FullFilePath=[FullFilePath,'shuffled\'] mkdir(FullFilePath) failedfilepath=[FullFilePath,'TooFewSpikes\']; mkdir(failedfilepath) K=0; order=1; fID=1; for h=1:1:length(NEVunits); if fID<=rf if NEVchannels(h)==Failed_timeframes{z,2}(fID,1)&NEVunits(h)==Failed_timeframes{z,2}(fID,2) fID=fID+1; s_SpikePhase=s_SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber); title([sprintf('(NA)(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(s_PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',s_PolarityProfile2(5,order+1))); hold off; saveas(K,[failedfilepath,sprintf('(NA)(Shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[failedfilepath,sprintf('(NA)(Shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end else s_SpikePhase=s_SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber); title([sprintf('(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(s_PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',s_PolarityProfile2(5,order+1))); hold off; saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end end else s_SpikePhase=s_SpikeFreqCell{3,h+1}; for j=1:1:length(LFPChannels); K=K+1; figure(K); subplot(2,1,1); H(K)=polarhistogram(s_SpikePhase(:,j),BinNumber); title([sprintf('(shuffled)%s %d-%ds Ch%dU%d LFP%d ',frametypetext,timeframe(1),timeframe(2),NEVchannels(h),NEVunits(h),LFPChannels(j)),Title1]); subplot(2,1,2); max_lim=1; x_fake=[0 max_lim 0 -max_lim]; y_fake=[max_lim 0 -max_lim 0]; h_fake=compass(x_fake,y_fake); hold on; compass(s_PolarityProfiles(5,order)); set(h_fake,'Visible','off'); title(sprintf('PI:%f',s_PolarityProfile2(5,order+1))); hold off; saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.fig',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); saveas(K,[FullFilePath,sprintf('(shuffled)SpikePhase_Chan%d_U%d_LFP%d_%d-%ds.jpg',NEVchannels(h),NEVunits(h),LFPChannels(j),timeframe(1),timeframe(2))]); order=order+1; end end end for t=1:rf A=find(str2double(s_PolarityProfile2(1,:))==Failed_timeframes{z,2}(t,1)); B=find(str2double(s_PolarityProfile2(2,A))==Failed_timeframes{z,2}(t,2)); s_PolarityProfile2(:,A(B))=[]; s_PolarityProfiles(:,A(B)-1)=[]; A=find(cell2mat(s_SpikeFreqCell(1,2:end))==Failed_timeframes{z,2}(t,1))+1; B=find(cell2mat(s_SpikeFreqCell(2,A))==Failed_timeframes{z,2}(t,2)); s_SpikeFreqCell(:,A(B))=[]; end filename=sprintf('(shuffled)SpikePhaseAnalysis_%d-%ds_%s.mat',timeframe(1),timeframe(2),frametypetext); save([FullFilePath,filename],'SourceFile','s_SpikeFreqCell','s_PolarityProfiles','s_PolarityProfile2'); compilation{2,z}=timeframe; compilation{3,z}=frametypetext; compilation{4,z}=PolarityProfile2; compilation{7,z}=s_PolarityProfile2; end if kilo==1 FullFilePath=[GetFilePath,'(kilo)SpikePhaseAnalysisTimeFrames(v2_3b6)\']; elseif kilo==0 FullFilePath=[GetFilePath,'SpikePhaseAnalysisTimeFrames(v2_3b6)\']; end save([FullFilePath,'PolarityCompilation.mat'],'SourceFile','compilation');