Browse Source

gin commit from med068pc183

Modified files: 9
Deleted files: 2
Andrey Vinogradov 2 years ago
parent
commit
a811a54d4f

+ 19 - 10
Codes/MATLAB spike detection/DetectSpikes_Amp.m

@@ -1,28 +1,37 @@
-function [index,  spikes, thr, thrmax, noise_std_detect, noise_std_sorted] = DetectSpikes_Amp(Data, handles)
+function [index,  spikes, thr, thrmax, noise_std_detect] = DetectSpikes_Amp(Data, handles)
 % DetectSpikes_Amp Detect spikes based on amplitude thresholding. Adapted
 % from WaveClus
-%   Detailed explanation goes here
-%   index - spike times in ms
-%   spikes - spike waveforms
+% INPUT
+%   Data - raw waveform from the electrode
+%   handles - detection parameters
+% OUTPUT
+%   index - spike times (the only output written into result .csv files)
+%   spikes - spike waveforms (not in use outside of the function currently)
+%   thr - minimum threshold used for detection (not in use outside of the function)
+%   thrmax - maximum threshold used for artefact removal (not in use outside of the function)
+%   noise_std_detect - estimate of the standard deviation of the background
+%   noise for detection (not in use outside of the function)
+
 
 index_all=[];
 spikes_all=[];
-% SPIKE DETECTION WITH AMPLITUDE THRESHOLDING
-% korjaa sampleratella
-% julkaisun defaultit: 24000Hz, pre 20 events, post 44 events
+
+%% ACTUAL SAMPLE RATE BASED PRE/POST correction and REFRACTORY PERIOD calculation
+% default: 24000Hz, pre 20 events, post 44 events, actual sr is read and
+% the values are corrected
 sr_factor = 24000;
 handles.w_pre = ceil(handles.w_pre/sr_factor*handles.sr);      %number of pre-event data points stored (def. 20)
 handles.w_post = ceil(handles.w_post/sr_factor*handles.sr);    %number of post-event data points stored (def. 44)
 handles.ref = floor(handles.min_ref_per*handles.sr/1000);      %number of counts corresponding to the dead time 
 
-
-[spikes,thr,thrmax, noise_std_detect, noise_std_sorted, index]  = amp_detect(Data,handles);       %detection with amp. thresh. --> voiko Data olla taulukkona?
+%% SPIKE DETECTION 
+[spikes,thr,thrmax, noise_std_detect, index]  = amp_detect(Data,handles);       %detection with amp. thresh. --> voiko Data olla taulukkona?
 
 index = index';
 
 index_all = [index_all index];
 spikes_all = [spikes_all; spikes];
-index = index_all/handles.sr;                  %spike times in ms --> now in s
+index = index_all/handles.sr;                 
 spikes = spikes_all;
 
 

+ 61 - 113
Codes/MATLAB spike detection/FileRunner.m

@@ -3,44 +3,39 @@ classdef FileRunner
     %   Detailed explanation goes here
     
     properties (SetAccess = private)
-        files; % lista tiedostoista, jotka luetaan
-        last_file; % viimeinen onnistunut, alotetaan 0:sta
-        file_read_error; % alotetaan 0:sta
-        
-       
-        
+        files; % listing of the files that are read
+        last_file; % the last file starts from 0
+        file_read_error; % reading error starts from 1     
         type; % type of file reading now
         spkdetconf; % configuration for detecting spikes
         recordtime; % timestamp for the current recording
      
-        duration; % duration of current file
-        
-        
+        duration; % duration of current file      
         samplerate;
         
-        file_finished; % alotetaan 1:st?
-        
-        % axionfileille OpenAxionFile-metodissa
-        first_ch_col; % alotus 1
-        final_ch_col; % luetaan maksimi tiedoston tiedoista
-        first_ch_row; % alotus 1
-        final_ch_row; % luetaan maksimi tiedoston tiedoista
-        first_well_col; % alotus 1
-        final_well_col; % luetaan maksimi tiedoston tiedoista
-        first_well_row; % alotus 1
-        final_well_row; % luetaan maksimi tiedoston tiedoista
-        last_ch_col; % alotus 0
-        last_ch_row; % alotus 0
-        last_well_row; % alotus 0
-        last_well_col; % alotus 0
-        
-        Plate_type
-        T; %added by andrey table for csv
+        file_finished; % finishing indicator starts from 1
+        
+        first_ch_col; % starts with 1
+        final_ch_col; % depends on plate type
+        first_ch_row; % starts with 1
+        final_ch_row; % depends on plate type
+        first_well_col; % starts with 1
+        final_well_col; % depends on plate type
+        first_well_row; % starts with 1
+        final_well_row; % depends on plate type
+        last_ch_col; % starts with 0
+        last_ch_row; % starts with 0
+        last_well_row; % starts with 0
+        last_well_col; % starts with 0
+        
+        Plate_type;
+        T; % table for csv
         csv_directory; %where to save the csv files
     end
     
     methods
         function obj = FileRunner(p_SpikeDetConfiguration, p_TargetFiles,csv_directory)
+            %constructor method
             obj.files = p_TargetFiles;
             obj.last_file = 0;
             
@@ -51,28 +46,14 @@ classdef FileRunner
             obj.csv_directory=csv_directory;
         end
         
-        function r_fname = getCurrentFile(obj)
-            % files; % lista tiedostoista, jotka luetaan
-            % last_file; % viimeinen onnistunut, alotetaan 0:sta
+        function r_fname = getCurrentFile(obj)    
             r_fname = obj.files{obj.last_file};
         end
-        
-        %function r_durr = getCurrentFileDuration(obj) %NOT IN USE ?!!!!!!!!!!!
-            %r_durr = obj.duration;
-        %end
-        
-        % onnistuiko,indexes,spikes,chlabel,duration,record_time,detection threshold,noise
-        function [r_result, obj] = ReadSpikesFromNextChannel(obj) % palauttaa seuraavan kanavan indexes & waves etc
-            %curr_well=strcat(char(64+obj.last_well_row+1),num2str(obj.last_well_col+1));
-            %splitter=@(x) strsplit(x, '/'); 
-            %while ~any(strcmp(cellfun(@(z) z{end}, cellfun(splitter,{h5info(obj.getCurrentFile()).Groups(1).Groups.Name}, 'UniformOutput',false), 'UniformOutput',false), curr_well)) 
-                %fprintf([curr_well '_' 'seems to be excluded\n']);
-                %obj = obj.JumpNextWell;
-                %curr_well=strcat(char(64+obj.last_well_row+1),num2str(obj.last_well_col+1));
-            %end
-            % onko viel? elektrodeita, jotka pit?isi lukea
-            if obj.file_finished &&...% jos oli viimeinen tai eka kanava
-                obj.last_file >= length(obj.files) % ja jos ei ole tiedostoja jotka pit?isi lukea                        
+
+        function [r_result, obj] = ReadSpikesFromNextChannel(obj) 
+            %proceeds the processing to the next channel/file
+            if obj.file_finished &&...
+                obj.last_file >= length(obj.files)                      
                 Finalcsv=sortrows(obj.T, 2);
                 thefilenamefull=strsplit(obj.getCurrentFile(), '.');
                 thefilename=strsplit(thefilenamefull{1}, '\');
@@ -80,8 +61,7 @@ classdef FileRunner
                 writetable(Finalcsv, strcat(obj.csv_directory, '\',thefilename,'_', 'spikes.csv'));
                 disp('All files have been read from!');
                 r_result = [];
-            else % ei ole eka tai vika kanava, tai on tiedostoja joista pit?isi lukea
-                % Tarvitseeko avata uusi tiedosto
+            else 
                 if obj.file_finished&&...
                     obj.last_file>0
                     Finalcsv=sortrows(obj.T, 2);
@@ -94,49 +74,33 @@ classdef FileRunner
                 end
                 
                 if obj.file_finished
-                    obj = obj.OpenNewFile(); %( obj.files(obj.last_file+1) );
+                    obj = obj.OpenNewFile(); 
                 end
                 
-                %if ~obj.file_read_error
-                %Reading the channel 
-
-                
-                [r_result.error,m_data, r_result.chlabel, obj] = obj.ReadAxionChannel();
-                
-                %disp(r_result.chlabel)
-                
+                [r_result.error,m_data, r_result.chlabel, obj] = obj.ReadAxionChannel();               
                 obj = obj.SiirraOsoittimetSeuraavaan();
                 if ~r_result.error
-                    % check if the channel is behaving properly or if is noisy  
                     obj.spkdetconf.sr = obj.samplerate; 
-                    [r_result.indexes, r_result.spikes, r_result.threshold, ~, r_result.noise, ~] = DetectSpikes_Amp(m_data, obj.spkdetconf);
-                                                           
+                    [r_result.indexes, r_result.spikes, r_result.threshold, ~, r_result.noise] = DetectSpikes_Amp(m_data, obj.spkdetconf);                                                           
+                    %only r_result.indexes are in use currently from all of the
+                    %output of DetectSpikes_Amp function
                     Time=r_result.indexes;                  
                     Channel=string(repmat(r_result.chlabel, numel(r_result.indexes),1));
                     %repeating the channel label as many times %as the number
                     %of the detected spikes. (For the table creation)
                     TT=table(Channel, Time);
                     obj.T=[obj.T; TT];
-                end
-                %end
+                end                
             end
         end
      
                
                  
         function [errorread,data,chlabel, obj] = ReadAxionChannel(obj)
-            
-            %well does not exist jump to another and display...might be not
-            %recorded or so
-            %jump to the next
-            %if channel (single)does not exist then display channel not recorded 
-            %next
+            %reading the data from the electrode
             disp(obj.getCurrentFile())
-            data = []; %there in Meeri's code to return at least something with e.g. error=1
-            %error = 0;
-            %duration = [];
-            m_well=strcat(char(64+obj.last_well_row+1),num2str(obj.last_well_col+1));
-            
+            data = [];
+            m_well=strcat(char(64+obj.last_well_row+1),num2str(obj.last_well_col+1));            
             splitter=@(x) strsplit(x, '/'); 
             while ~any(strcmp(cellfun(@(z) z{end}, cellfun(splitter,{h5info(obj.getCurrentFile()).Groups(1).Groups.Name}, 'UniformOutput',false), 'UniformOutput',false), m_well)) 
                 fprintf([m_well '_' 'seems to be excluded\n']);
@@ -150,36 +114,30 @@ classdef FileRunner
             end
             
             m_ch=strcat(num2str(obj.last_ch_col+1),num2str(obj.last_ch_row+1));
-            chlabel = strcat(m_well,'_',m_ch); % which well which electrode
+            chlabel = strcat(m_well,'_',m_ch); % form the label (well and electrode)
             disp(['now reading: ' chlabel]);
-            try 
-                
+            try                 
                 data = h5read(obj.getCurrentFile(),['/Data/' m_well '/' m_ch]);
-                %if isempty(obj.duration)
-                    %duration = length(data)/obj.samplerate;
-                    %obj.duration = duration;
-                %else
-                    %duration = obj.duration;%better when open file?????
-                %end
                 errorread = 0;
-            catch %err
+            catch 
                 errorread = 1;
                 InactiveChannels=h5read(obj.getCurrentFile(), '/DataInfo/InactiveChannels');
                 if any(contains(InactiveChannels, chlabel)) 
                     disp(['Channel' ' ' chlabel ' ' 'is inactive, it is not recorded in the initial analysed file.'])
                 else
-                    error(['Smth wrong with' ' ' chlabel ' ' 'channel reading!']); %ADDED
+                    error(['Smth wrong with' ' ' chlabel ' ' 'channel reading!']); 
                 end
           
             end 
         end
         
         function obj = SiirraOsoittimetSeuraavaan(obj)
-            if obj.final_ch_row-1 <= obj.last_ch_row % viimeinen elektrodi, otetaan seuraavalta sarakkeelta
+            %shifts the well/electrode addresses (row/col) to the next after processing 
+            if obj.final_ch_row-1 <= obj.last_ch_row 
                 obj.last_ch_row = obj.first_ch_row-1;
-                if obj.final_ch_col-1 <= obj.last_ch_col % viimeinen sarake, otetaan seuraava kaivo
+                if obj.final_ch_col-1 <= obj.last_ch_col 
                     obj.last_ch_col = obj.first_ch_col-1;
-                    if obj.final_well_col-1 <= obj.last_well_col % viimeinen kaivosarake, otetaan seuraavalta rivilt?
+                    if obj.final_well_col-1 <= obj.last_well_col 
                         obj.last_well_col = obj.first_well_col-1;
                         if obj.final_well_row-1 <= obj.last_well_row
                             disp('All channels have been read from!');
@@ -194,40 +152,31 @@ classdef FileRunner
                     obj.last_ch_col = obj.last_ch_col+1;
                 end                    
             else
-                obj.last_ch_row = obj.last_ch_row+1; % riveitt?in seuraava electrodi
-            end
-            
+                obj.last_ch_row = obj.last_ch_row+1; 
+            end            
         end
         
         
         function obj = OpenNewFile(obj)
+            %launches new file opening
             obj.file_finished = 0;
-            obj.file_read_error = 1; %setting it to 1 if everything successful it will be switched to 0 during OpenAxionFile
-            while obj.file_read_error && obj.last_file < length(obj.files) %&& obj.spkdetconf.MinTime > obj.duration
+            obj.file_read_error = 1; %setting it to 1, if everything successful it will be switched to 0 during OpenAxionFile
+            while obj.file_read_error && obj.last_file < length(obj.files) 
                 disp(['Opening new file: ' obj.files{obj.last_file+1}]);
                 obj.type = strsplit(obj.files{obj.last_file+1},'.');
                 obj.type = cell2mat(obj.type(end));
-                                                 
-                %strcmp('raw',obj.type) % axion files
-                obj = obj.OpenAxionFile(obj.files{obj.last_file+1});
-              
-                obj.last_file = obj.last_file+1;
-                
+                obj = obj.OpenAxionFile(obj.files{obj.last_file+1});            
+                obj.last_file = obj.last_file+1;                
             end
             if obj.file_read_error && obj.last_file == length(obj.files)
                 error('NOT A SINGLE FILE WAS OPENED SUCCESSFULLY!')
             end
         end
-        
-
-  
-        
-        function obj = OpenAxionFile(obj, p_TargetFile) % sets also sample rate, voltageScale, recordtime
+                 
+        function obj = OpenAxionFile(obj, p_TargetFile) 
+            % sets the sample rate, voltageScale, recordtime from the file
             try
-                %obj.AxisFileData=AxisFile(p_TargetFile);
-               
-                obj.Plate_type=h5readatt(p_TargetFile, '/DataInfo/', 'Plate type');
-                
+                obj.Plate_type=h5readatt(p_TargetFile, '/DataInfo/', 'Plate type');                
                 obj.first_ch_col = 1;
                 obj.first_ch_row = 1;
                 obj.first_well_col = 1;
@@ -253,18 +202,17 @@ classdef FileRunner
                 end
                  
                 obj.duration = h5readatt(p_TargetFile, '/DataInfo/','DurationInSec');
-                obj.samplerate = h5readatt(p_TargetFile, '/DataInfo/','SamplingFrequencyInHz');
-                %obj.recordtime = ...;
-                
+                obj.samplerate = h5readatt(p_TargetFile, '/DataInfo/','SamplingFrequencyInHz');               
                 obj.file_read_error = 0;
-            catch %err
+            catch 
                 disp (['File ' p_TargetFile ' could not be opened successfully!']);
                 obj.file_read_error = 1;
             end
         end
         
         
-        function obj = JumpNextWell(obj)        
+        function obj = JumpNextWell(obj)  
+            %proceeds to the next well 
             if obj.final_well_col-1 <= obj.last_well_col 
                 obj.last_well_col = obj.first_well_col-1;
                 if obj.final_well_row-1 <= obj.last_well_row

+ 6 - 9
Codes/MATLAB spike detection/Main.m

@@ -20,22 +20,19 @@ parameters.w_post = 44;                    %number of post-event data points sto
 % pre + post ==> amount of points stored per spike
 
 parameters.detection = 'both';             %type of threshold 'both','neg' or 'pos'
-parameters.stdmin = 4.5;                   %minimum threshold (def. 5) for the operation with SWTTEO 4.5 should be selected
-parameters.stdmax = 50;                    %maximum threshold. Is used for artefact removal
-parameters.interpolation = 'y';            %interpolation for alignment
-parameters.int_factor = 2;                 %interpolation factor (def. 2)
-parameters.detect_fmin = 200;              %high pass filter for detection (def. 300)
+parameters.stdmin = 4.5;                   %minimum threshold factor (def. 5) for the operation with SWTTEO 4.5 should be selected
+parameters.stdmax = 50;                    %maximum threshold factor. Is used for artefact removal
+parameters.interpolation = 'y';            %interpolation for alignment (not in use here as spike storage is not involved)
+parameters.int_factor = 2;                 %interpolation factor (interpolation for spike storing is not in use)
+parameters.detect_fmin = 200;              %high pass filter for detection (def. 200)
 parameters.detect_fmax = 3000;             %low pass filter for detection (def. 3000)
-parameters.sort_fmin = 200;                %high pass filter for sorting (def. 300)
-parameters.sort_fmax = 3000;               %low pass filter for sorting (def. 3000)
 parameters.segments = 1;                   %nr. of segments in which the data is cutted. --> segmentit voisi analysoida rinnakkain?
-%handles.par.sr = [];                        %sampling frequency in Hz extracted from the recording file metadata during the processiong
 
 %setting of the output directory 
 csv_directory = uigetdir('title','SELECT YOUR OUTPUT FOLDER');
 %% RUN
 runrunman = FileRunner(parameters, list(:,1),csv_directory); %constructor method application
-[result, runrunman] = runrunman.ReadSpikesFromNextChannel(); 
+[result, runrunman] = runrunman.ReadSpikesFromNextChannel(); %reading the next channel
 
 while ~isempty(result)
         

+ 31 - 60
Codes/MATLAB spike detection/amp_detect.m

@@ -1,15 +1,14 @@
-function [spikes,thr,thrmax, noise_std_detect, noise_std_sorted, index] = amp_detect(x,handles)
+function [spikes,thr,thrmax, noise_std_detect, index] = amp_detect(x,handles)
 % Detect spikes with amplitude thresholding. Uses median estimation.
 % Detection is done with filters set by fmin_detect and fmax_detect. Spikes
 % are stored for sorting using fmin_sort and fmax_sort. This trick can
 % eliminate noise in the detection but keeps the spikes shapes for sorting.
 
-%The initial (only amp thresholding) version was taken by Meeri from
-%some ?? Wave Clus PROGRAM Get_spikes ?? or so
+%The initial version involves only amp thresholding 
 %Stationary Wavelet Transform Teager Energy Operator block is added by Andrey.
 %Workflow: %1. Amplitude thresholding detects spikes in the prefiltered signal with TCF=4.5
            %2. The number of the detected spikes is fed into SWTTEO algorithm
-           %3. The spike detected by both algorithms are considered as True Positives (within error window)        
+           %3. The spike detected by both algorithms are considered as True Positives (within an error window)        
                
 x = double(x);
 
@@ -18,9 +17,6 @@ x = double(x);
 % other               :  filtfilt
 % other               :  fix_filter
 % other               :  int_spikes
-%
-% 090814 korjattu ohi-indeksointi lis??m?ll? tarkistukset spike storing
-% osioon
 
 sr=handles.sr;
 w_pre=handles.w_pre;
@@ -31,41 +27,27 @@ stdmin = handles.stdmin;
 stdmax = handles.stdmax;
 fmin_detect = handles.detect_fmin;
 fmax_detect = handles.detect_fmax;
-fmin_sort = handles.sort_fmin;
-fmax_sort = handles.sort_fmax;
 index1 = [];
 spikes = [];
 thr = [];
 thrmax = [];
 noise_std_detect = [];
-noise_std_sorted = [];
 
-% HIGH-PASS FILTER OF THE DATA
 
-xf=zeros(length(x),1);
-if exist('ellip')                         %Checks for the signal processing toolbox
-    [b,a]=ellip(2,0.1,40,[fmin_detect fmax_detect]*2/sr);
-    xf_detect=filtfilt(b,a,x);
-    [b,a]=ellip(2,0.1,40,[fmin_sort fmax_sort]*2/sr);
-    xf=filtfilt(b,a,x);
-else
-    xf=fix_filter(x);                   %Does a bandpass filtering between [300 3000] without the toolbox.
-    xf_detect = xf;
-end
-lx=length(xf);
+%% HIGH-PASS FILTER OF THE DATA
+%Checks for the signal processing toolbox
+[b,a]=ellip(2,0.1,40,[fmin_detect fmax_detect]*2/sr);
+xf=filtfilt(b,a,x);
 
-%clear x; MOVED TO THE END 
+noise_std_detect = median(abs(xf))/0.6745;
 
-noise_std_detect = median(abs(xf_detect))/0.6745;
-noise_std_sorted = median(abs(xf))/0.6745;
-thr = stdmin * noise_std_detect;        %thr for detection is based on detected settings.
-thrmax = stdmax * noise_std_sorted;     %thrmax for artifact removal is based on sorted settings.
-
-% LOCATE SPIKE TIMES
+thr = stdmin * noise_std_detect;        %thr for detection
+thrmax = stdmax * noise_std_detect;     %thrmax for artifact removal 
+%% LOCATE SPIKE TIMES
 switch detect
     case 'pos'
         nspk = 0;
-        xaux = find(xf_detect(w_pre+2:end-w_post-2) > thr) +w_pre+1;
+        xaux = find(xf(w_pre+2:end-w_post-2) > thr) +w_pre+1;
         xaux0 = 0;
         for i=1:length(xaux)
             if xaux(i) >= xaux0 + ref
@@ -77,7 +59,7 @@ switch detect
         end
     case 'neg'
         nspk = 0;
-        xaux = find(xf_detect(w_pre+2:end-w_post-2) < -thr) +w_pre+1;
+        xaux = find(xf(w_pre+2:end-w_post-2) < -thr) +w_pre+1;
         xaux0 = 0;
         for i=1:length(xaux)
             if xaux(i) >= xaux0 + ref
@@ -89,10 +71,7 @@ switch detect
         end
     case 'both'
         nspk = 0;
-        % palauttaa ei-nollat ja lis?? niihin w_pre+1
-        % absoluuttisista arvoista
-        % jotka xf_detectiss? v?lill? w_pre+2:end-w_post-2
-        xaux = find(abs(xf_detect(w_pre+2:end-w_post-2)) > thr) +w_pre+1;
+        xaux = find(abs(xf(w_pre+2:end-w_post-2)) > thr) +w_pre+1;
         xaux0 = 0;
         for i=1:length(xaux)
             if xaux(i) >= xaux0 + ref
@@ -101,18 +80,19 @@ switch detect
                 if isempty(iaux)
                     iaux = 0;
                 end
-                index1(nspk) = iaux + xaux(i) -1; % heitt?? virheilmon: index sis?lt?? spiketimet, maxi ja iaux tyhji?
+                index1(nspk) = iaux + xaux(i) -1; 
                 xaux0 = index1(nspk);
             end
         end
 end
-%some values for SWTTEO
-in.M = xf_detect;
+%% some values for SWTTEO
+in.M = xf;
 in.SaRa = sr;
 %Parameters
-params.method = 'numspikes';
-params.numspikes = numel(index1);
-params.filter = 0;
+params.method = 'numspikes'; %detection with predefined number of spikes
+params.numspikes = numel(index1); %same number of spikes as detected with the amplitude thresholding
+params.filter = 0; %no addition filtering for SWTTEO algorithm
+%Wavelet decomposition level selection based on the sampling frequency
 if sr==10000
     params.wavLevel=2;
 elseif sr==12500
@@ -125,11 +105,10 @@ else
     error('Is it smth else than MCS or Axion? I know only sampling frequencies 10000, 12500, 25000 and 50000!')
 end
 
-%Getting result for the SET NUMSPIKES from SWTTEO algorithm
+%% Getting result for the SET NUMSPIKES from SWTTEO algorithm
 index2= SWTTEO(in,params);
 
-%Selecting the spikes detected by both algorithms within some []error
-%interval
+%% Selecting the spikes detected by both algorithms within some []error interval
 errortime=ceil(2e-4*sr);%(in samples)
 
 check=zeros(1, params.numspikes);
@@ -142,29 +121,20 @@ index=index(index>0);
 index=sort(index);
  
 nspk=numel(index);
-% SPIKE STORING (with or without interpolation)
-% tehd??n vain jos spikej? indexiss?!
+%% SPIKE STORING (with or without interpolation) 
+% [SPIKE STORING OUTPUT IS NOT IN USE CURRENTLY, ONLY ARTIFACT REMOVAL]  
+
 if ~isempty(index)
     ls=w_pre+w_post;
     spikes=zeros(nspk,ls+4);
-    xf=[xf zeros(1,w_post)];              % xf ja zeros ei yhteensopivia
+    xf=[xf zeros(1,w_post)];              
     parfor i=1:nspk                          %Eliminates artifacts
         if max(abs( xf(index(i)-w_pre:index(i)+w_post) )) < thrmax
-            % jos piikin indeksi datan vika tai tokavika piste niin aiheuttaa
-            % ohi indeksoinnin, eik? piikkimuotokaan ole kokonainen
-            % eli eliminoidaan artifactina
             if ~(index(i)+w_post+2 > length(xf))
-                % lis?t??n nollaborderia, koska piikki aivan datan alussa
                 if index(i)-w_pre-1 < 1 
                     spikes(i,:)=[0 xf(abs(index(i)-w_pre:index(i)+w_post+2))];
                 else
                     spikes(i,:)=xf(abs(index(i)-w_pre-1:index(i)+w_post+2));
-                    % aiheuttaa ongelman jos index(i)-w_pre-1 < 1
-                    % tai
-                    % jos index(i)+w_post+2 > length(xf)
-                    % eli jos piikki on datassa alle w_pren p??ss? alusta
-                    % tai aivan datan lopussa
-                    % --> t?ll?in ei saada kokonaista piikkimuotoa, joten
                 end
             end
         end
@@ -173,15 +143,16 @@ if ~isempty(index)
     spikes(aux,:)=[];
     index(aux)=[];
     
-    clear x; %moved from line 57
+    clear x; 
 
 end
         
 switch handles.interpolation
+    %spike storing is currently not in use
     case 'n'
         %eliminates borders that were introduced for interpolation 
-        spikes(:,end-1:end)=[];       % laittaa kaikille vikan ja tokavikan []
-        spikes(:,1:2)=[];             % laittaa kaikille ekan ja tokan []
+        spikes(:,end-1:end)=[];       
+        spikes(:,1:2)=[];             
     case 'y'
         %Does interpolation
         spikes = int_spikes(spikes,handles);   

+ 0 - 34
Codes/MATLAB spike detection/fix_filter.m

@@ -1,34 +0,0 @@
-function y = fix_filter(x)
-% fix_filter (in case signal processing toolbox is not available).
-% filters data x with an eliptic passband between [300 3000] Hz.
-
-a = [1.0000 -2.3930  2.0859 -0.9413 0.2502];
-b = [0.1966 -0.0167 -0.3598 -0.0167 0.1966];
-
-x = x(:);  
-len = size(x,1);  
-b = b(:).';
-a = a(:).';
-nb = length(b);
-na = length(a);
-nfilt = max(nb,na);
-
-nfact = 3*(nfilt-1);  % length of edge transients
-
-rows = [1:nfilt-1  2:nfilt-1  1:nfilt-2];
-cols = [ones(1,nfilt-1) 2:nfilt-1  2:nfilt-1];
-data = [1+a(2) a(3:nfilt) ones(1,nfilt-2)  -ones(1,nfilt-2)];
-sp = sparse(rows,cols,data);
-zi = sp \ ( b(2:nfilt).' - a(2:nfilt).'*b(1) );
-
-y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];
-
-y = filter(b,a,y,[zi*y(1)]);
-y = y(length(y):-1:1);
-y = filter(b,a,y,[zi*y(1)]);
-y = y(length(y):-1:1);
-
-y([1:nfact len+nfact+(1:nfact)]) = [];
-
-y = y.';   
-

+ 1 - 0
Codes/MATLAB spike detection/int_spikes.m

@@ -1,5 +1,6 @@
 function [spikes1] = int_spikes(spikes,handles); 
 %Interpolates with cubic splines to improve alignment.
+%Not in use currently 
 
 w_pre=handles.w_pre;
 w_post=handles.w_post;

+ 8 - 5
Codes/PCA_code/PCA.m

@@ -1,23 +1,26 @@
-%PCA analysis script written by Laura Ylä-Outinen
+%% PCA analysis script written by Laura Ylä-Outinen
+
+%Data organisation
+
 X1=readtable('PCA_table.xlsx');
-X1(1:36,1)=table({'hPSC_20517_MEA1and2'});
+X1(1:36,1)=table({'hPSC_20517_MEA1and2'}); %merging MEA1 and MEA2 data as belonging to the same experiment
 X2=X1(:,4:end);
 X3=table2array(X2);
 Xstandard=zscore(X3);
 expressions=Xstandard;
-%%
+%% subtipes selection
 
 subtypes=X1(:,1);
 subtypes= table2array(subtypes);
 subtypes3=categorical(subtypes);
 
 
-%%
+%% PCA calculation
 [coeff, score, latent, tsquared, explained,mu] = pca(expressions); %if matrix do not contain blanks
 X=score(:,1);
 Y=score(:,2);
 Z=score(:,3);
 
-%%
+%% 3D Visualisation with "gscatter3b" function
 gscatter3b(X,Y,Z, subtypes3, [0.82 0.3 0.1; 0.82 0.3 0.1; 0.82 0.3 0.1; 0 0.48 0.83; 0 0.48 0.83; 0 0.48 0.83;] ,'ox.ox.',[10,10,15,10,10,15]) ;
 

+ 0 - 1
Codes/meaRtools/.Rhistory

@@ -1 +0,0 @@
-source('D:/gin-cli-latest-windows64/attempt2/Codes/meaRtools/MEA_analysis_Axion.R')

+ 2 - 4
Codes/meaRtools/MEA_analysis_Axion.R

@@ -2,7 +2,7 @@
 # Analysis code for MEA data
 
 #library(IGM.MEA)
-library(meaRtools) #attach meaRtools after IGM.MEA to get right parameter names <- fixed by getting rid of IGM.MEA, just by creating an extra function
+library(meaRtools) #attach meaRtools after IGM.MEA to get right parameter names <- fixed by getting rid of IGM.MEA, just by creating an extra .R function
 library(plyr)
 library(ggplot2)
 library(reshape2)
@@ -12,8 +12,6 @@ library(pracma)
 library(svDialogs) 
 
 
-CRAPels=list()
-
 #CODEdirectory<- folder where the code files reside
 CODEdirectory<-paste0(choose.dir(default = "", caption = "SELECT THE FOLDER WITH THE CODE"),"\\")
 #alternative, not to select this every code launch:
@@ -80,7 +78,7 @@ source(paste0(CODEdirectory, "github_nb_bursts_modifications_LAa.R"))
 
 source(paste0(CODEdirectory, "cheminfo.R"))
 
-
+CRAPels=list()
 ########## IMPORT ##########
 
 my.output <- paste0(output_folder, '/Analysis')

+ 1 - 0
Codes/meaRtools/cheminfo.R

@@ -1,4 +1,5 @@
 cheminfo<-function (file, masterChemFile = masterChemFile) 
+  #function from IGM.MEA package 
 {
   masterChem = read.csv(masterChemFile)
   masterCD <- as.data.frame(masterChem)

+ 19 - 1
README.md

@@ -83,7 +83,25 @@ Codes directory contains the following folders:
 * Connectivity Analysis
 * PCA code 
 
-The implementational details can be found in the original publication.
+
+
+## Code Usage Notes
+To launch the MATLAB spike detection code, the user needs to open the *Main.m* file and then add the whole folder containing the MATLAB analysis code to the MATLAB path (so that the program is capable of finding the functions that the code requires). Next, by pressing MATLAB’s green “Run” button, the analysis is launched. The selection of the .h5 raw data files for the analysis is implemented via a pop-up window that opens as the user launches the code. One or more files can be selected for analysis. The output folder selection is performed in the same way. The code sequentially analyses the selected files, providing the user with spike .csv files as an output. If the user decides to use different parameters than described in the methods section for spike detection, the parameters to be changed are located in *Main.m* and *amp_detect.m*.
+To implement the analyses in meaRtools, the user needs to follow the steps specified below:
+i. First, the user needs to open *MEA_analysis_Axion.R* file.
+ii. Then, by clicking the “Source” button, the code is launched. The user sees pop-up windows for selecting the code-containing folder, the output folder, the spike .csv files, the noisy electrode file and the expLog file.
+iii. The last pop-up window enables the selection of the analyzed MEA type. There are 12- and 48-well MEA plates available, and the user only needs to enter an integer that corresponds to the analyzed plate type.
+The folder selection windows sometimes do not appear on top of the RStudio window; then, they are found in the Windows taskbar. If the user wants to avoid the code directory selection step, it is possible to remove the first pop-up window by replacing the first "choose.dir" function with the code directory address.
+It should be noted that electrodes that
+* have no detected spikes (are not mentioned in the MATLAB-generated .csv spike files)
+* are listed in the noisy electrode .csv file
+* are eliminated by the minimum-spikes-per-minute criterion
+are not included in the meaRtools analysis.
+Essentially, if for these abovementioned reasons all electrodes in a particular well are cancelled for all DIVs included in the analysis, this well is not displayed in the analysis output files.
+There is a possibility of implementing the code in segments by performing segment selection and pressing the “Run” button.
+The PCA and connectivity analysis MATLAB code packages are delivered in their corresponding folders. To run the PCA code, the user needs to open the *PCA.m* code in the "PCA" folder and add the whole folder to the MATLAB path. The folder contains a table with preselected activity features of the cell populations obtained during the provided analysis path. The next step is to click the "Run" button to launch PCA.
+To launch the connectivity analysis, the user needs to open the *Connectivity.m* file in the "Connectivity analysis" folder and add this folder to the MATLAB path. After clicking the "Run" button, the pop-up window for .h5 file selection appears. After selecting the desired file, a new pop-up window for MEA well selection appears. The user needs to take into account the list of excluded wells, which are automatically displayed in the command window after the file selection, and avoid selecting them. For the connectivity analysis, the threshold value for connectivity strength can be changed in the script *cross_selection_correlated_channels.m*. More information on the CorSE and analysis guidelines can be found at https://se.mathworks.com/matlabcentral/fileexchange/59626-spectral-entropy-based-neuronal-network-synchronization-analysis-corse.
+
 
 ## Related Publications
 * Kapucu, F.E., Vinogradov A., Hyvärinen T., Ylä-Outinen L., Narkilahti S. Comparative microelectrode array dataset of functional development of human PSC-derived and rat neuronal networks. Submitted