ConvertPhy2KlustaSuite does not suit phy files

[eliezyer]

This function doesn't work with the new phy files (spike_clusters.npy and cluster_group.tsv). It needs an update

[brendonw1]

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley notifications@github.com wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#event-1483459195, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

[petersenpeter]

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson notifications@github.com:

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley notifications@github.com wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#event-1483459195, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367213108, or mute the thread https://github.com/notifications/unsubscribe-auth/ANb8uVsvKoHvDDTdJno_UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc .

[brendonw1]

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen notifications@github.com wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson notifications@github.com:

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley <notifications@github.com

wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#event-1483459195, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe-auth/ANb8uVsvKoHvDDTdJno_UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367322635, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc .

[eliezyer]

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson notifications@github.com wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen <notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson notifications@github.com:

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#event-1483459195, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe-auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367322635 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367323602, or mute the thread https://github.com/notifications/unsubscribe-auth/AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml.
% % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels')
%if ~isempty(par.SpkGrps(groupidx).Channels)
% for each group loop through, find all templates clus
tgroup          = allgroups(groupidx);%shank number
ttemplateidxs   = find(templateshankassignments==tgroup);%which templates/clusters are in that shank
ttemplates      = templates(:,:,ttemplateidxs);
tPCFeatureInds  = pcFeatureInds(:,ttemplateidxs);

tidx            = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank
tclu            = clu;%extract template/cluster assignments of spikes on this shank
tspktimes       = spktimes;

gidx            = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group
channellist     = [];

temp.SpkGrps(1).Channels = 0:63;
for ch = 1:length(temp.SpkGrps)
    if ismember(gidx(1),temp.SpkGrps(ch).Channels+1)
        channellist = temp.SpkGrps(ch).Channels+1;
        break
    end
end
if isempty(channellist)
    disp(['Cannot find spkgroup for group ' num2str(groupidx) ])
    continue
end

%% spike extraction from dat
if groupidx == 1;
    dat             = memmapfile(datpath,'Format','int16');
end
tsampsperwave   = (sbefore+safter);
ngroupchans     = length(channellist);
valsperwave     = tsampsperwave * ngroupchans;
wvforms_all     = zeros(length(tspktimes)*tsampsperwave*ngroupchans,1,'int16');
wvranges        = zeros(length(tspktimes),ngroupchans);
wvpowers        = zeros(1,length(tspktimes));

for j=1:length(tspktimes)
    try
        w = dat.data((double(tspktimes(j))-sbefore).*totalch+1:(double(tspktimes(j))+safter).*totalch);
        wvforms=reshape(w,totalch,[]);
        %select needed channels
        wvforms = wvforms(channellist,:);
%         % detrend
%         wvforms = floor(detrend(double(wvforms)));
        % median subtract
        wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter);
        wvforms = wvforms(:);

    catch
        disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']);
        disp(['Time range of that spike was: ' num2str(double(tspktimes(j))-sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples'])
        wvforms = zeros(valsperwave,1);
    end

    %some processing for fet file
    wvaswv = reshape(wvforms,tsampsperwave,ngroupchans);
    wvranges(j,:) = range(wvaswv);
    wvpowers(j) = sum(sum(wvaswv.^2));

    lastpoint = tsampsperwave*ngroupchans*(j-1);
    wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms;
%     wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)')));
    if rem(j,100000) == 0
        disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done'])
    end
end
wvranges = wvranges';

%% Spike features

% for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file
fets    = zeros(sum(tidx),size(pcFeatures,2),ngroupchans);
pct     = repmat(pcFeatures,1,1,5);
pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group
%and rearrange their features to match the shank order
allshankclu = unique(tclu);

for tc = 1:length(allshankclu)
    tsc     = allshankclu(tc);
    i       = find(tclu==tsc);
    tforig  = pct(i,:,:);%the subset of spikes with this clu ide
    tfnew   = tforig; %will overwrite

    ii      = tdx(tc,:);%handling nan cases where the template channel used was not in the shank
    gixs    = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels
    bixs    = isnan(ii);
    g       = ii(gixs);

    tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements
    tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank
    try
        fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels));
    catch
        keyboard
    end
end
%extract for relevant spikes only...
% and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless
tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx).Channels)));%lazy reshaping
tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx).Channels)));
tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx).Channels)));
fets = cat(2,tfet1,tfet2,tfet3)';%     fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']);

% fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(tspktimes')); fets = cat(1,double(fets),double(wvpowers),double(wvranges),double(tspktimes')); fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]);
resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]);
fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]);
spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]);

%fet SaveFetIn(fetname,fets);

%clu
% if ~exist(cluname,'file')
    tclu = cat(1,length(unique(tclu)),double(tclu));
    fid=fopen(cluname,'w'); 
%     fprintf(fid,'%d\n',clu);
    fprintf(fid,'%.0f\n',tclu);
    fclose(fid);
    clear fid
% end

%res
fid=fopen(resname,'w'); 
fprintf(fid,'%.0f\n',tspktimes);
fclose(fid);
clear fid

%spk
fid=fopen(spkname,'w'); 
fwrite(fid,wvforms_all,'int16');
fclose(fid);
clear fid 

disp(['Shank ' num2str(tgroup) ' done'])
%end
%end

end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf 
    BufInd = [(i-1)*nBuf+1:min(i*nBuf,size(Fet,1))];
    temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))];
    fprintf(outputfile,formatstring,temp');
end

end fclose(outputfile);

[brendonw1]

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer notifications@github.com wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson notifications@github.com wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson notifications@github.com:

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe-auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367322635 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,'int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))-sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx).Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx).Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx).Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges),double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367365604, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

[petersenpeter]

The problem arises as the original list of peak channels (calculated by Kilosort) is not updated when new units are created in phy. So I just did a plugin for phy which exports the shank number for each unit to a .npy file. I will update the phy plugins with these changes.

2018-02-21 10:43 GMT-05:00 Brendon Watson notifications@github.com:

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer notifications@github.com wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson < notifications@github.com> wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson <notifications@github.com :

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 https://github.com/buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe- auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#

issuecomment-367322635

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,' int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))- sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx). Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx).Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx).Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/ ' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers, wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges), double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367365604 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367368406, or mute the thread https://github.com/notifications/unsubscribe-auth/ANb8uXDYET6wtNb4ydY6jTOKE6kdDxfmks5tXDmLgaJpZM4SMuyc .

[brendonw1]

Fabulous. Thank you

On Wed, Feb 21, 2018 at 4:04 PM, Peter C Petersen notifications@github.com wrote:

The problem arises as the original list of peak channels (calculated by Kilosort) is not updated when new units are created in phy. So I just did a plugin for phy which exports the shank number for each unit to a .npy file. I will update the phy plugins with these changes.

2018-02-21 10:43 GMT-05:00 Brendon Watson notifications@github.com:

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer notifications@github.com wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson < notifications@github.com> wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson < notifications@github.com :

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 https://github.com/ buzsakilab/buzcode/issues/149 to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe- auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#

issuecomment-367322635

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);% find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,' int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))- sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx). Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx).Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx).Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/ ' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers, wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges), double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367365604 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367368406 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ANb8uXDYET6wtNb4ydY6jTOKE6kdDxfmks5tXDmLgaJpZM4SMuyc .

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367471655, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTdtZkoqyPwkTphRwAnS9Jpq9mbAtks5tXITpgaJpZM4SMuyc .

[brendonw1]

Let me know if I should help. On Wed, Feb 21, 2018 at 4:06 PM Brendon Watson brendon.watson@gmail.com wrote:

Fabulous. Thank you

On Wed, Feb 21, 2018 at 4:04 PM, Peter C Petersen < notifications@github.com> wrote:

The problem arises as the original list of peak channels (calculated by Kilosort) is not updated when new units are created in phy. So I just did a plugin for phy which exports the shank number for each unit to a .npy file. I will update the phy plugins with these changes.

2018-02-21 10:43 GMT-05:00 Brendon Watson notifications@github.com:

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer notifications@github.com wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson < notifications@github.com> wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson < notifications@github.com :

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 < https://github.com/buzsakilab/buzcode/issues/149> to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe- auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#

issuecomment-367322635

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,' int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))- sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx). Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx).Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx).Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/ ' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers, wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges), double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367365604 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub < https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367368406>, or mute the thread < https://github.com/notifications/unsubscribe-auth/ANb8uXDYET6wtNb4ydY6jTOKE6kdDxfmks5tXDmLgaJpZM4SMuyc

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367471655, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTdtZkoqyPwkTphRwAnS9Jpq9mbAtks5tXITpgaJpZM4SMuyc .

[eliezyer]

I've made the script attached to convert the phy files to neurosuite. Can you guys try and see if everything is ok?

On Wed, Feb 21, 2018 at 9:28 PM, Brendon Watson notifications@github.com wrote:

Let me know if I should help. On Wed, Feb 21, 2018 at 4:06 PM Brendon Watson brendon.watson@gmail.com wrote:

Fabulous. Thank you

On Wed, Feb 21, 2018 at 4:04 PM, Peter C Petersen < notifications@github.com> wrote:

The problem arises as the original list of peak channels (calculated by Kilosort) is not updated when new units are created in phy. So I just did a plugin for phy which exports the shank number for each unit to a .npy file. I will update the phy plugins with these changes.

2018-02-21 10:43 GMT-05:00 Brendon Watson notifications@github.com:

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer notifications@github.com wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson < notifications@github.com> wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson < notifications@github.com :

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 < https://github.com/buzsakilab/buzcode/issues/149> to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe- auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#

issuecomment-367322635

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,' int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))- sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx). Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx). Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx). Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/ ' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers, wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges), double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367365604 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub < https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367368406 , or mute the thread < https://github.com/notifications/unsubscribe-auth/ ANb8uXDYET6wtNb4ydY6jTOKE6kdDxfmks5tXDmLgaJpZM4SMuyc

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367471655, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTdtZkoqyPwkTphRwAnS9Jpq9mbAtks5tXITpgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367546735, or mute the thread https://github.com/notifications/unsubscribe-auth/AIKHOASAICplyJiWHnmDoptEZHXUlmJWks5tXNDqgaJpZM4SMuyc .

function ConvertPhyKilo2Neurosuite(basepath,basename)

% Converts Phy kilosort output files into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat, .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % It assumes that your directory contains: % spike_clusters.npy, spikes_times.npy, cluster_group.tsv, pcfeatures.npy, % templates.npy, cluster_ids.npy and shanks.npy (these last two are % generated by a phy plugin -Export Shanks- found in the lab wiki) % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name) % % Eliezyer de Oliveira, 2018

savepath = basepath;

if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% getting spike information

spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); % pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))'; templates = readNPY('templates.npy');

mkdir(fullfile(savepath,'PhyClus')) %% assigning cluster ids to shanks shanks = readNPY('shanks.npy'); cluShank = readNPY('cluster_ids.npy');

auxC = unique(clu); templateshankassignments = zeros(size(auxC)); for idx = 1:length(auxC) temp = find(cluShank == auxC(idx)); templateshankassignments(idx) = shanks(temp); end grouplookup = rez.ops.kcoords; allgroups = unique(grouplookup); allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels')
%if ~isempty(par.SpkGrps(groupidx).Channels)
% for each group loop through, find all templates clus
tgroup          = allgroups(groupidx);%shank number
ttemplateidxs   = find(templateshankassignments==tgroup);%which templates/clusters are in that shank

% ttemplates = templates(:,:,ttemplateidxs); % tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx            = ismember(clu,auxC(ttemplateidxs));%find spikes indices in this shank
tclu            = clu(tidx);%extract template/cluster assignments of spikes on this shank
tspktimes       = spktimes(tidx);

gidx            = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group
channellist     = [];

for ch = 1:length(par.SpkGrps)
    if ismember(gidx(1),par.SpkGrps(ch).Channels+1)
        channellist = par.SpkGrps(ch).Channels+1;
        break
    end
end
if isempty(channellist)
    disp(['Cannot find spkgroup for group ' num2str(groupidx) ])
    continue
end

%% spike extraction from dat
if groupidx == 1;
    dat             = memmapfile(datpath,'Format','int16');
end
tsampsperwave   = (sbefore+safter);
ngroupchans     = length(channellist);
valsperwave     = tsampsperwave * ngroupchans;
wvforms_all     = zeros(length(tspktimes)*tsampsperwave*ngroupchans,1,'int16');
wvranges        = zeros(length(tspktimes),ngroupchans);
wvpowers        = zeros(1,length(tspktimes));

for j=1:length(tspktimes)
    try
        w = dat.data((double(tspktimes(j))-sbefore).*totalch+1:(double(tspktimes(j))+safter).*totalch);
        wvforms=reshape(w,totalch,[]);
        %select needed channels
        wvforms = wvforms(channellist,:);
%         % detrend
%         wvforms = floor(detrend(double(wvforms)));
        % median subtract
        wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter);
        wvforms = wvforms(:);

    catch
        disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']);
        disp(['Time range of that spike was: ' num2str(double(tspktimes(j))-sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples'])
        wvforms = zeros(valsperwave,1);
    end

    %some processing for fet file
    wvaswv = reshape(wvforms,tsampsperwave,ngroupchans);
    wvranges(j,:) = range(wvaswv);
    wvpowers(j) = sum(sum(wvaswv.^2));

    lastpoint = tsampsperwave*ngroupchans*(j-1);
    wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms;
%     wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)')));
    if rem(j,100000) == 0
        disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done'])
    end
end
wvranges = wvranges';

%% Spike features

% for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:length(ttemplateidxs); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file
fets    = zeros(sum(tidx),size(pcFeatures,2),ngroupchans);
pct     = pcFeatures(tidx,:,:);

%for each cluster/template id, grab at once all spikes in that group
%and rearrange their features to match the shank order
allshankclu = unique(tclu);

for tc = 1:length(allshankclu)
    tsc     = allshankclu(tc);
    i       = find(tclu==tsc);
    tforig  = pct(i,:,:);%the subset of spikes with this clu ide
    tfnew   = tforig; %will overwrite

    ii      = tdx(tc,:);%handling nan cases where the template channel used was not in the shank
    gixs    = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels
    bixs    = isnan(ii);
    g       = ii(gixs);

    tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements
    tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank
    try
        fets(i,:,:) = tfnew(:,:,1:size(pct,3));
    catch
        keyboard
    end
end
%extract for relevant spikes only...
% and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless
tfet1 = squeeze(fets(:,1,1:size(pct,3)));%lazy reshaping
tfet2 = squeeze(fets(:,2,1:size(pct,3)));
tfet3 = squeeze(fets(:,3,1:size(pct,3)));
fets = cat(2,tfet1,tfet2,tfet3)';%     fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']);

% fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(tspktimes')); fets = cat(1,double(fets),double(wvpowers),double(wvranges),double(tspktimes')); fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]);
resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]);
fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]);
spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]);

%fet SaveFetIn(fetname,fets);

%clu
% if ~exist(cluname,'file')
    tclu = cat(1,length(unique(tclu)),double(tclu));
    fid=fopen(cluname,'w'); 
%     fprintf(fid,'%d\n',clu);
    fprintf(fid,'%.0f\n',tclu);
    fclose(fid);
    clear fid
% end

%res
fid=fopen(resname,'w'); 
fprintf(fid,'%.0f\n',tspktimes);
fclose(fid);
clear fid

%spk
fid=fopen(spkname,'w'); 
fwrite(fid,wvforms_all,'int16');
fclose(fid);
clear fid 

disp(['Shank ' num2str(tgroup) ' done'])
%end
%end

end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'PhyClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf 
    BufInd = [(i-1)*nBuf+1:min(i*nBuf,size(Fet,1))];
    temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))];
    fprintf(outputfile,formatstring,temp');
end

end fclose(outputfile);

[brendonw1]

Sorry where can we find it? On Fri, Feb 23, 2018 at 6:22 PM eliezyer notifications@github.com wrote:

I've made the script attached to convert the phy files to neurosuite. Can you guys try and see if everything is ok?

On Wed, Feb 21, 2018 at 9:28 PM, Brendon Watson notifications@github.com wrote:

Let me know if I should help. On Wed, Feb 21, 2018 at 4:06 PM Brendon Watson <brendon.watson@gmail.com

wrote:

Fabulous. Thank you

On Wed, Feb 21, 2018 at 4:04 PM, Peter C Petersen < notifications@github.com> wrote:

The problem arises as the original list of peak channels (calculated by Kilosort) is not updated when new units are created in phy. So I just did a plugin for phy which exports the shank number for each unit to a .npy file. I will update the phy plugins with these changes.

2018-02-21 10:43 GMT-05:00 Brendon Watson notifications@github.com:

We ended up doing something like ook at the max channel and figuring out which shank that was on and then assigning it to that shank... does that work?

On Wed, Feb 21, 2018 at 10:34 AM, eliezyer < notifications@github.com> wrote:

So the temporary solution for me have been to put all the cluster in only one shank using script attached (it's based on ConvertKilosort2Neurosuite). I'm using spike_times.npy, spike_cluster.npy and cluter_group.tsv. I can't get the shank/channel info from templates because they are from original kilosort clusters and are not the same after spike sorting on phy.

On Wed, Feb 21, 2018 at 8:21 AM, Brendon Watson < notifications@github.com> wrote:

Can you send me your script and I’ll see if I can do it?

Also maybe we can pull in the .xml or metadata files if necessary but hopefully we don’t have to use that. On Wed, Feb 21, 2018 at 8:17 AM Peter C Petersen < notifications@github.com

wrote:

I have my own script for loading the units into Matlab. Eli was trying to do a script to export the phy output to clu and rez files, but he is having trouble with getting the shanks info out.

2018-02-20 23:52 GMT-05:00 Brendon Watson < notifications@github.com :

Right. I actually don't know the format of the new phy files. Does anyone know? Have Peter or anyone done this? Peter is using the new phy, does he know how to read out the files? If you can give some pointers I can try to write the code

On Tue, Feb 20, 2018 at 4:56 PM, David Tingley < notifications@github.com

wrote:

Assigned #149 < https://github.com/buzsakilab/buzcode/issues/149> to @brendonw1 https://github.com/brendonw1.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# event-1483459195 , or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZEYXRdPh9Ba8UxfHtl6guZuspLJks5tWz-HgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367213108 , or mute the thread < https://github.com/notifications/unsubscribe- auth/ANb8uVsvKoHvDDTdJno_ UUlDTvrPO-aBks5tW6EYgaJpZM4SMuyc

.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#

issuecomment-367322635

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTWnLGRZT1splAINjl5imNnq-YjtUks5tXBdwgaJpZM4SMuyc

.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367323602 , or mute the thread https://github.com/notifications/unsubscribe-auth/ AIKHONrYOMSbs0cTfy5JYoGodM03Cil6ks5tXBg-gaJpZM4SMuyc .

function ConvertNPY2Neurosuite_KSW(rez)

% Converts KiloSort templates Klusta into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat and an .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name)

% Brendon Watson 2016 % Edited by Peter Petersen 2017

savepath = cd; basepath = cd; basename = rez.ops.basename; if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% % [spikeTimes, ii] = sort(spikeTimes); spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); amplitudes = readNPY('amplitudes.npy'); amplitudes = amplitudes(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';

mkdir(fullfile(savepath,'SortedClus')) %% do homework for assigning templates to shanks % [~,shank]=fileparts(basepath); templates = permute(readNPY('templates.npy'),[3 2 1]); % m = min(templates,[],2);%find the min value of each waveform on each channel % [~,m] = min(m,[],1);%find which channel minimum is least % m = squeeze(m);%which channel is minimum on each template m = max(abs(templates),[],2);%find the most deviated value of each waveform on each channel [~,m] = max(m,[],1);%find which channel has most deviated value for each templnate m = squeeze(m);%squeeze to 1d vector

grouplookup = rez.ops.kcoords;%list of group/shank of each channel templateshankassignments = grouplookup(m);%for the list of maximal channels, which group is each in allgroups = unique(grouplookup); allgroups = 1; %Grp 0 contain discared channels allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank ttemplates = templates(:,:,ttemplateidxs); tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(auxiliarC,auxiliarC);%ismember(clu,ttemplateidxs);%find spikes indices in this shank tclu = clu;%extract template/cluster assignments of spikes on this shank tspktimes = spktimes;

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

temp.SpkGrps(1).Channels = 0:63; for ch = 1:length(temp.SpkGrps) if ismember(gidx(1),temp.SpkGrps(ch).Channels+1) channellist = temp.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,' int16');

wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double( tspktimes(j))+safter).totalch); wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))- sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:size(ttemplates,3); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = repmat(pcFeatures,1,1,5); pct = cat(3,pct,pct(:,:,1:4));

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:length(temp.SpkGrps(groupidx).Channels)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:length(temp.SpkGrps(groupidx). Channels)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:length(temp.SpkGrps(groupidx). Channels))); tfet3 = squeeze(fets(:,3,1:length(temp.SpkGrps(groupidx). Channels))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/ ' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers, wvranges,double(tspktimes'));

fets = cat(1,double(fets),double(wvpowers),double(wvranges), double(tspktimes'));

fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'SortedClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <

https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367365604

, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTZ1B5qIR4F8MOBXdyhC8hIlMLEzqks5tXDeXgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <

https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367368406 ,

or mute the thread < https://github.com/notifications/unsubscribe-auth/ ANb8uXDYET6wtNb4ydY6jTOKE6kdDxfmks5tXDmLgaJpZM4SMuyc

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149# issuecomment-367471655, or mute the thread https://github.com/notifications/unsubscribe-auth/ ADXrTdtZkoqyPwkTphRwAnS9Jpq9mbAtks5tXITpgaJpZM4SMuyc .

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <https://github.com/buzsakilab/buzcode/issues/149#issuecomment-367546735 , or mute the thread < https://github.com/notifications/unsubscribe-auth/AIKHOASAICplyJiWHnmDoptEZHXUlmJWks5tXNDqgaJpZM4SMuyc

.

function ConvertPhyKilo2Neurosuite(basepath,basename)

% Converts Phy kilosort output files into klusters-compatible % fet,res,clu,spk files. Works on a single shank of a recording, assumes a % 16bit .dat, .xml file is present in "basepath" (home folder) and % that they are named basename.dat and basename.xml. % It assumes that your directory contains: % spike_clusters.npy, spikes_times.npy, cluster_group.tsv, pcfeatures.npy, % templates.npy, cluster_ids.npy and shanks.npy (these last two are % generated by a phy plugin -Export Shanks- found in the lab wiki) % % Inputs: % basepath - directory path to the main recording folder with .dat and .xml % as well as shank folders made by makeProbeMapKlusta2.m (default is % current directory matlab is pointed to) % basename - shared file name of .dat and .xml (default is last part of % current directory path, ie most immediate folder name) % % Eliezyer de Oliveira, 2018

savepath = basepath;

if ~exist('basepath','var') [~,basename] = fileparts(cd); basepath = cd; end if ~exist('rez','var'); load(fullfile(basepath,'rez.mat')) end

Nchan = rez.ops.Nchan; connected = rez.connected; xcoords = rez.xc; ycoords = rez.yc; % Nchan = rez.ops.Nchan; % connected = ones(Nchan, 1); % xcoords = ones(Nchan, 1); % ycoords = (1:Nchan)';

par = LoadXml(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels; sbefore = 16;%samples before/after for spike extraction safter = 24;%... could read from SpkGroups in xml if isfield(par.SpkGrps,'nSamples') if ~isempty(par.SpkGrps(1).nSamples); if isfield(par.SpkGrps,'PeakSample') if ~isempty(par.SpkGrps(1).PeakSample); sbefore = par.SpkGrps(1).PeakSample; safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample; end end end end

if exist(rez.ops.fbinary,'file') datpath = rez.ops.fbinary; else datpath = fullfile(basepath,[basename '.dat']); end

%% identify timestamps that are from good clusters spikeTimes = readNPY('spike_times.npy'); clusters = readNPY('spike_clusters.npy'); S = tdfread('cluster_group.tsv'); group = S.group; cluster_id = S.cluster_id; c = 0; for idx = 1:size(group,1) if group(idx,1:4) == 'good' c = c+1; GClusters(c) = idx; end end ExtClus = cluster_id(GClusters);

% Separating idx by cluster auxiliarC = []; for idx = 1:length(ExtClus) auxC = find(clusters == ExtClus(idx)); auxiliarC = cat(1,auxiliarC,auxC); end [~,I] = sort(auxiliarC); auxiliarC = auxiliarC(I); %% getting spike information

spktimes = uint64(readNPY('spike_times.npy')); spktimes = spktimes(auxiliarC); clu = uint32(readNPY('spike_clusters.npy')); clu = clu(auxiliarC); pcFeatures = readNPY('pc_features.npy'); pcFeatures = pcFeatures(auxiliarC,:,:); % pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))'; templates = readNPY('templates.npy');

mkdir(fullfile(savepath,'PhyClus')) %% assigning cluster ids to shanks shanks = readNPY('shanks.npy'); cluShank = readNPY('cluster_ids.npy');

auxC = unique(clu); templateshankassignments = zeros(size(auxC)); for idx = 1:length(auxC) temp = find(cluShank == auxC(idx)); templateshankassignments(idx) = shanks(temp); end grouplookup = rez.ops.kcoords; allgroups = unique(grouplookup); allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

%if isfield(par.SpkGrps(groupidx),'Channels') %if ~isempty(par.SpkGrps(groupidx).Channels) % for each group loop through, find all templates clus tgroup = allgroups(groupidx);%shank number ttemplateidxs = find(templateshankassignments==tgroup);%which templates/clusters are in that shank % ttemplates = templates(:,:,ttemplateidxs); % tPCFeatureInds = pcFeatureInds(:,ttemplateidxs);

tidx = ismember(clu,auxC(ttemplateidxs));%find spikes indices in this shank tclu = clu(tidx);%extract template/cluster assignments of spikes on this shank tspktimes = spktimes(tidx);

gidx = 1:length(rez.ops.kcoords);%find(rez.ops.kcoords == tgroup);%find all channels in this group channellist = [];

for ch = 1:length(par.SpkGrps) if ismember(gidx(1),par.SpkGrps(ch).Channels+1) channellist = par.SpkGrps(ch).Channels+1; break end end if isempty(channellist) disp(['Cannot find spkgroup for group ' num2str(groupidx) ]) continue end

%% spike extraction from dat if groupidx == 1; dat = memmapfile(datpath,'Format','int16'); end tsampsperwave = (sbefore+safter); ngroupchans = length(channellist); valsperwave = tsampsperwave ngroupchans; wvforms_all = zeros(length(tspktimes)tsampsperwave*ngroupchans,1,'int16'); wvranges = zeros(length(tspktimes),ngroupchans); wvpowers = zeros(1,length(tspktimes));

for j=1:length(tspktimes) try w = dat.data((double(tspktimes(j))-sbefore).totalch+1:(double(tspktimes(j))+safter).totalch);

wvforms=reshape(w,totalch,[]); %select needed channels wvforms = wvforms(channellist,:); % % detrend % wvforms = floor(detrend(double(wvforms))); % median subtract wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter); wvforms = wvforms(:);

catch disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']); disp(['Time range of that spike was: ' num2str(double(tspktimes(j))-sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples']) wvforms = zeros(valsperwave,1); end

%some processing for fet file wvaswv = reshape(wvforms,tsampsperwave,ngroupchans); wvranges(j,:) = range(wvaswv); wvpowers(j) = sum(sum(wvaswv.^2));

lastpoint = tsampsperwavengroupchans(j-1); wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms; % wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)'))); if rem(j,100000) == 0 disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done']) end end wvranges = wvranges';

%% Spike features % for each template, rearrange the channels to reflect the shank order tdx = []; for tn = 1:length(ttemplateidxs); % tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template for k = 1:length(channellist); % i = find(tTempPCOrder==channellist(k)); if ~isempty(k) tdx(tn,k) = k; else tdx(tn,k) = nan; end end end

featuresperspike = 3; % kilosort default

% initialize fet file fets = zeros(sum(tidx),size(pcFeatures,2),ngroupchans); pct = pcFeatures(tidx,:,:);

%for each cluster/template id, grab at once all spikes in that group %and rearrange their features to match the shank order allshankclu = unique(tclu);

for tc = 1:length(allshankclu) tsc = allshankclu(tc); i = find(tclu==tsc); tforig = pct(i,:,:);%the subset of spikes with this clu ide tfnew = tforig; %will overwrite

ii = tdx(tc,:);%handling nan cases where the template channel used was not in the shank gixs = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels bixs = isnan(ii); g = ii(gixs);

tfnew(:,:,gixs) = tforig(:,:,g);%replace ok elements tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank try fets(i,:,:) = tfnew(:,:,1:size(pct,3)); catch keyboard end end %extract for relevant spikes only... % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless tfet1 = squeeze(fets(:,1,1:size(pct,3)));%lazy reshaping tfet2 = squeeze(fets(:,2,1:size(pct,3))); tfet3 = squeeze(fets(:,3,1:size(pct,3))); fets = cat(2,tfet1,tfet2,tfet3)';% fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']); % fets = double(squeeze(fets(1,:,:))); %mean activity per spike % fetmeans = mean(fets,1);%this is pretty redundant with wvpowers % %find first pcs, make means of those... % featuresperspike = 4; % firstpcslist = 1:featuresperspike:size(fets,1); % firstpcmeans = mean(fets(firstpcslist,:),1); % % nfets = size(fets,1)+1; % fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(tspktimes')); fets = cat(1,double(fets),double(wvpowers),double(wvranges),double(tspktimes')); fets = fets'; % fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk

cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]); resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]); fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]); spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]); %fet SaveFetIn(fetname,fets);

%clu % if ~exist(cluname,'file') tclu = cat(1,length(unique(tclu)),double(tclu)); fid=fopen(cluname,'w'); % fprintf(fid,'%d\n',clu); fprintf(fid,'%.0f\n',tclu); fclose(fid); clear fid % end

%res fid=fopen(resname,'w'); fprintf(fid,'%.0f\n',tspktimes); fclose(fid); clear fid

%spk fid=fopen(spkname,'w'); fwrite(fid,wvforms_all,'int16'); fclose(fid); clear fid

disp(['Shank ' num2str(tgroup) ' done']) %end %end end clear dat copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'PhyClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize) BufSize = inf; end

nFeatures = size(Fet, 2); formatstring = '%d'; for ii=2:nFeatures formatstring = [formatstring,'\t%d']; end formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w'); fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))]; fprintf(outputfile,formatstring,temp'); else nBuf = floor(size(Fet,1)/BufSize)+1;

for i=1:nBuf BufInd = [(i-1)nBuf+1:min(inBuf,size(Fet,1))]; temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))]; fprintf(outputfile,formatstring,temp'); end end fclose(outputfile);

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/buzsakilab/buzcode/issues/149#issuecomment-368167564, or mute the thread https://github.com/notifications/unsubscribe-auth/ADXrTbm08SPTC5KgRp0qzKaV1ytSMlG-ks5tX0gigaJpZM4SMuyc .

[eliezyer]

In this link: https://github.com/eliezyer/buzcode/blob/master/utilities/ConvertPhyKilo2Neurosuite.m

[DavidTingley]

Once you guys are happy with this, there is a PR (#151 ) that will incorporate Elie's code into buzcode/dev

[lisaroux]

Hi all, thanks Eli for writing this script. Has anyone tested it since then?

[DavidTingley]

Closing this, conversation moved to #225