Browse Source

gin commit from med068pc183

New files: 1
Modified files: 6
Andrey Vinogradov 2 years ago
parent
commit
7a889fc3b5

+ 8 - 11
Codes/Connectivity Analysis/Connectivity.m

@@ -28,6 +28,7 @@
 
 
 TargetFile=[path file];
+%Checking plate type and selecting well/channel characteristics accordingly
 Plate_type=h5readatt(TargetFile, '/DataInfo/', 'Plate type');
 if strcmp(Plate_type, '12-well plate')
     final_ch_col = 8;
@@ -41,28 +42,25 @@ elseif strcmp(Plate_type, '48-well plate')
     final_well_row = 6;
 else
     disp('Unknown MEA layout!')
-    
-    
-    
+   
 end
 
+%% Reading ExcludedWells & InactiveChannels
 ExcludedWells=h5read(TargetFile, '/DataInfo/ExcludedWells');
 disp('Note the excluded wells!:')
 disp(ExcludedWells)
 InactiveChannels=h5read(TargetFile, '/DataInfo/InactiveChannels');
 disp('Inactive channels!:')
 disp(InactiveChannels)
-%%
+%% Well selection
 prompt='Select your well.Please, use capital letters, e.g. A1:';
 MyWell = inputdlg(prompt);
 MyWell=MyWell{:};
-%%
-well_row=double(MyWell(1)-64);
 
+well_row=double(MyWell(1)-64);
 well_col=str2double(MyWell(2:end));
 
-
-%%
+%% Checking if the well data is present there in the file
 splitter=@(x) strsplit(x, '/'); 
 RecordedWells={h5info(TargetFile).Groups(1).Groups.Name};
 if ~any(strcmp(cellfun(@(z) z{end}, cellfun(splitter,RecordedWells, 'UniformOutput',false), 'UniformOutput',false), MyWell)) 
@@ -70,8 +68,7 @@ if ~any(strcmp(cellfun(@(z) z{end}, cellfun(splitter,RecordedWells, 'UniformOutp
 end
  
  
- %%
- 
+ %% Reading channels in the selected well
   
 for i = 1 : final_ch_row
       
@@ -89,5 +86,5 @@ for i = 1 : final_ch_row
     end
 end
    
-
+%% Moving the data forward for processing. Launching "Spectral_entropy_HighPass_data_adaptive.m" function 
     Spectral_entropy_HighPass_data_adaptive;

+ 5 - 5
Codes/Connectivity Analysis/OnlyCross_data_adaptive_zero_delay.m

@@ -1,10 +1,9 @@
-%%CROSS-CORRELATION
+%% CROSS-CORRELATION
 cross=zeros(chNO,chNO);
 delay=zeros(chNO,chNO);
-CorrLength = (2*size(A,2))-1
+CorrLength = (2*size(A,2))-1;
 meancorr=zeros(chNO,CorrLength);
 
-
 % ALL THE CHANNELS
 for channel=1:chNO
     
@@ -25,7 +24,7 @@ for channel=1:chNO
         
         [xcf,lags,bounds]=crosscorr(y1,y2,(size(A,2)-1));
         
-        Correlations(i,1:CorrLength)=xcf;%calculating cross-correlations' mean for each channel
+        Correlations(i,1:CorrLength)=xcf;%calculating cross-correlations for each channel
         
         maxDelayforCorr = 0;% max delay for max correlation in seconds
         zeroLagindex = find(lags==0);
@@ -63,5 +62,6 @@ for channel=1:chNO
     
 end
 
-%%
+%% Moving forward for the threshold setting and selection of the sufficient correlations
+% "with cross_selection_correlated_channels.m" function
 cross_selection_correlated_channels()

+ 27 - 43
Codes/Connectivity Analysis/Spectral_entropy_HighPass_data_adaptive.m

@@ -1,5 +1,5 @@
 %% ENTROPY ALGORITHM --> LOADING ONE CHANNEL
-fs = 12500; %sampling freq.
+fs = 12500; %sampling frequency (for Axion)
 WS = fs*0.5; %window size to calculate Spectral Entropy in samples fs for 1s
 overlap=0.5;%other than 0.5 makes it more complex and code needs to be changed
 TotalDataSample = length(DataCell{1,7});%length data in samples
@@ -13,62 +13,46 @@ for channel=1:chNO
     channel
     if ~isempty(DataCell{channel,7})
         x=DataCell{channel,7}(1:end);%Travel MEA has problems with the first index otherwise start from 1
-%     x=x-mean(x); %If you subtract mean from entire signal
-    %% LOADING SIGNAL WINDOWS EQUAL TO 1 SECOND
-    
-        for i=1:SElength-1 
-        %% THEORETICAL SPECTRAL RESOLUTION OF 1 Hz-OVERLAP = 50%
-        
+        %     x=x-mean(x); %If you subtract mean from entire signal
+        %% LOADING SIGNAL WINDOWS EQUAL TO 1 SECOND
+        for i=1:SElength-1
+            %% THEORETICAL SPECTRAL RESOLUTION OF 1 Hz-OVERLAP = 50%
             y=x((i-1)*WS*overlap +1:(i-1)*WS*overlap+WS);
-        
             y=y-mean(y);%subtract mean from each window "preferable"
-      
-        %%
-        %**************************************************************
-        % GET VECTORS b AND a HERE (use function rico.m)
-        %**************************************************************
+            %%
+            %**************************************************************
+            % GET VECTORS b AND a HERE (use function rico.m)
+            %**************************************************************
             [b,a]=rico(0.01,6,50,1/fs);  % rico(z,B,fc,T)
-        %	z  (0.01) attenuazione minima alla frequenza fc
-        %   minimum attenuation at the center frequency frequency band
-
-        %	B  (2) larghezza di banda corrispondente alla attenuazione 0.707
-        %	Bandwidth which corresponds to attenuation 0.707
-
-        % 	fc (50) frequenza di centro banda 
-        %   Central frequency
-
-        % 	T  intervallo di campionamento 
-        %   Sampling interval 
-                     
+            %input parameters
+            %	z  (0.01) minimum attenuation at the center frequency frequency band
+            %	B  (2) Bandwidth which corresponds to attenuation 0.707
+            % 	fc (50) Central frequency
+            % 	T Sampling interval
+            % output parameters
+            %	b  filter coefficients (MA part)
+            %	a  filter coefficients (AR part)
             
             z=filtfilt(b,a,y);          % IIR
-        
-        % HIGH FILTER
+            %% HIGH FILTER
             n=3;
             wn=7/(fs/2); %7HzHPF
             [b,a]=butter(n,wn,'high');
-     
-        
             z1=filtfilt(b,a,z);
-        
-        %% POWER SPECTRUM OF FILTERED SIGNAL WITH RICURSIVE NOTCH (IIR) --> NO TIME DELAY
-        
+            
+            %% POWER SPECTRUM OF FILTERED SIGNAL WITH RICURSIVE NOTCH (IIR) --> NO TIME DELAY
             nfft=2^15;
             w=hann(WS);
             o=0;
-            [Pxx,f]=pwelch(z1,w,o,nfft,fs); 
-        
-            pxx_norm1=Pxx/sum(Pxx);%
-        
-        
-        
-        %% SPECTRAL ENTROPY (EQ.5)
-        
+            [Pxx,f]=pwelch(z1,w,o,nfft,fs);
+            pxx_norm1=Pxx/sum(Pxx);
+            
+            %% SPECTRAL ENTROPY (EQ.5)
             k=0;
             s1(i)=0;
             n=length(pxx_norm1);
             for j=1:n
-                if isfinite(log10(1/pxx_norm1(j)))==1 
+                if isfinite(log10(1/pxx_norm1(j)))==1
                     s1(i)=s1(i)+ pxx_norm1(j)*log10(1/pxx_norm1(j));
                 else
                     k=1;
@@ -76,7 +60,7 @@ for channel=1:chNO
                 end
             end
             s1(i)=s1(i)/log10(n);
-        
+            
         end
         A(channel,:)=s1;
     else
@@ -86,7 +70,7 @@ for channel=1:chNO
     
     
 end
-%%
+%%  Moving further for cross-correlation calculation with "OnlyCross_data_adaptive_zero_delay.m" function
 OnlyCross_data_adaptive_zero_delay()
 
 

+ 9 - 16
Codes/Connectivity Analysis/cross_selection_correlated_channels.m

@@ -1,3 +1,4 @@
+%Threshold setting and selection of the sufficient correlations for the visualisation 
 clear row
 clear col
 clear source
@@ -8,19 +9,16 @@ for t = 1 : length(cross)
 end
 %% eliminate NaN values
 cross(find(isnan(cross)))=0;
-%%
 
-% %% REMOVE ARTEFACTED CHANNELS
-% artefacted_channels = [4 19 37 46 62 63];
+%% REMOVE ARTEFACTED CHANNELS
+% artefacted_channels = [];
 % for u = 1 : length(artefacted_channels)
 %     cross(artefacted_channels(u) , :)=0;
 %     cross(: , artefacted_channels(u))=0;
 % end
-% %% 
-
-
-%SELECT TRESHOLD VALUE%
+%% SELECT TRESHOLD VALUE FOR CONNECTIVITY STRENGTH
 
+%HERE THE USER CAN CHANGE THE THRESHOLD VALUE FOR CONNECTIVITY STRENGTH
 treshold=0.7; %for fixed threshold
 [row,col]=find(cross>treshold); %for fixed threshold
 
@@ -28,13 +26,8 @@ if isempty(row)
     error('No correlations were found with the defined threshold value! Change the threshold value in cross_selection_correlated_channels.m')
 end
 
-
-
-
 j=0;
-for i=1:length(row)
-    
-    
+for i=1:length(row)   
     %         if delay(row(i),col(i)) >= 0;
     j=j+1;
     %             row(i)
@@ -42,12 +35,12 @@ for i=1:length(row)
     %             pause
     formatSpec = 'channels: %2.0f and %2.0f \n';
     fprintf(formatSpec,col(i),row(i))
-    source(j)=col(i);
-    target(j)=row(i);
+    source(j)=col(i); %source channels
+    target(j)=row(i); %target channels
     weight(j)=cross(row(i),col(i));
     %         end
     
 end
 N=length(source);
-%%
+%% Proceeding towards visualisation using "map.m" function
 map()

+ 10 - 19
Codes/Connectivity Analysis/map.m

@@ -1,4 +1,6 @@
-%% for one way connections comment out if not used
+%% VISUALISATION SCRIPT FOR CONNECTIVITY ANALYSIS
+
+%preparing source matrix, target matrix and weight matrix
 source_1=source;
 target_1=target;
 for i=1:length(source)
@@ -14,7 +16,6 @@ end
 source=(nonzeros(source))';
 target=(nonzeros(target))';
 weight=(nonzeros(weight))';
-%%
 
 clear channels
 clear ids
@@ -40,37 +41,28 @@ weight_new=zeros(N,N);
 
 
 for i=1:length(channels)
-    [row,col]=find(source==channels(i)) %make connections from source to target
+    [row,col]=find(source==channels(i)); %make connections from source to target
     if size(row) ~= 0
         x=target(row,col);
         x=x(1,:);
         for j=1: length(row)
-            [row1,col1]=find(channels==x(j))
+            [row1,col1]=find(channels==x(j));
             cm(i,col1)=1;
-            weight_new (i,col1) = cross(channels(i),channels(col1));%%THIS PART IS ADDED AFTER REALIZING THE MISTAKE IN PUBLICATION...REORDERS THE WEIGHTS OF LINKS RIGHT WAY
+            weight_new (i,col1) = cross(channels(i),channels(col1));%THIS REORDERS THE WEIGHTS OF LINKS RIGHT WAY
         end
     end
-   
 end
 
-
-
-
 for k=1:N
     ChannelName{k} = DataCell{channels(k),1}(end-3:end);
     ids{k}=['' ChannelName{k} ''''];
 end
 
-
-
-%%
+%% BIOGRAPH creates a bioinformatics graph object
 bg2=biograph(cm,ids,'NodeAutoSize','off','ShowArrows','off');
-
-    get(bg2.nodes,'ID')
-    set(bg2, 'EdgeType', 'curved')
-    hbg= view(bg2);
-
-
+get(bg2.nodes,'ID')
+set(bg2, 'EdgeType', 'curved')
+hbg= view(bg2);
 
 %% DO LAYOUT ACCORDING TO MEA SETUP
 page_Size = [400 400];
@@ -102,7 +94,6 @@ for t = 1 : length(weight_new)
        weight_col=[weight_col b];
    end
 end
-% added the correct till here
 
 for i = 1:length(weight_row)
 set(hbg.Edges(i),'LineWidth', (10*(abs(weight_new(weight_row(i),weight_col(i)))-treshold)+0.1)); % default (10*(abs(weight_new(weight_row(i),weight_col(i)))-treshold)+0.1 ))

+ 5 - 13
Codes/Connectivity Analysis/rico.m

@@ -1,19 +1,11 @@
-% filtro ricorsivo reiezione 50 Hz
 %recursive band rejection filter 50 Hz
-function [cMA,cAR]=rico(z,B,fc,T);
+function [cMA,cAR]=rico(z,B,fc,T)
 
 % input parameters
-%	z  (0.01) attenuazione minima alla frequenza fc
-%   minimum attenuation at the center frequency frequency band
-
-%	B  (2) larghezza di banda corrispondente alla attenuazione 0.707
-%	Bandwidth which corresponds to attenuation 0.707
-
-% 	fc (50) frequenza di centro banda 
-%   Central frequency
-
-% 	T  intervallo di campionamento 
-%   Sampling interval 
+%	z  (0.01) minimum attenuation at the center frequency frequency band fc
+%	B  (2) Bandwidth which corresponds to attenuation 0.707
+% 	fc (50) Central frequency
+% 	T Sampling interval 
 
 % output parameters
 %	cMA  filter coefficients (MA part)

+ 1 - 0
Codes/meaRtools/.Rhistory

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