Cluster Analysis and Quantification of Segmented Volumes from IMARIS 3D Segmentation Data

    % the goal of this script is to determine cluster information starting from
    % the results extracted from IMARIS 3D segmentation 
    clear all
    clast_par=1.3; 
    % filename = '20241107 1457 18h_0002-1.tif DAPI.xls';
    % filename = '202406 1457dsbp pRB473_1.xls';
    % filename = '1457 1h 2.xls';
    %filename = '20231214 1457dsbp pRB473 sbp alfa 508_003-1.xls';
    filename = '20240321 1457dsbppRB473sbp ALFA508 0 1.xls';
    %filename = '202406 1457dsbp pRB473_1.xls';
    %The only two input necessary are the name of the file to be analized i.e file name and parameter that control how far away two
    % segmented volumes should be. clast_par < 1 means the distance between two segmented volume should be inferior to the sum of their mean radius
    % clast_par =1 the distance between two segmented volumes should be equal to the sum of their mean radius
    %  clast_par >1 the distance between two segmented volumes should be bigger than the sum of their mean radius
    % the data to be quantified is stored in an excel multipage file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% segmented volumes radius determination
    %read_volume=xlsread(filename,32);% opens page where the volumes are listed and extracts the volume information for each cluster
    read_volume=readtable(filename, 'Sheet', 32);
    siz_point=size(read_volume);     %it extracts the information about how many volumes were quantified 
    radius_nn(1:siz_point(1))=0;     %it creates a vector that will contain the estimated radius of all the detected volumes
    %radius_nn=1000*( read_volume(:,1)*3/(4*pi)).^(1/3);
    radius_nn = 1000 * (read_volume.Volume * 3 / (4 * pi)).^(1 / 3);
    % the radius is calculated approximating the quantified volumes with be the volumes of a sphere
    Mean_radius_nm = mean(radius_nn)   % displays the mean radius of the segmented volumes
    %please note that the radius of the segmented volumes is a good
    %approximation of the full width at half maximum of the fluorescence object
    Variance_radius_n = std(radius_nn)  % displays the variance of the radius of the segmented volumes
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% determination of
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the distance between the segmented volumes and sorting them between volumes belonging to the same cluster and volumes not belonging to the same cluster.
    % center_of_mass=xlsread(filename,5);% open page where the coordinate of the center of mass of the segmented volumes are listed
    % 
    % total_num=size(center_of_mass); % again number of detected spots just to check if it match with the previous values (not necessary)
    % total_resolved_spots = total_num(1) % displays the total number of volumes quantified
    % %center_of_mass=xlsread(filename,5); % open the page of the file where the values of the centre of mass of the detected volumes are stored and write them into a matrix
    % distances(1:siz_point(1),1:siz_point(1))=0;% produce the matrix were the distances between detected volumes will be stored
    % cluster(1:siz_point(1),1:siz_point(1))=0;%produce a matrix where will be stored the information wether two segmented volumes are close together (distance fulfil cluster parameter) cluster (i,j)=1
    % %or not cluster(i,j)=0;
    % 
    % %%calculates the matrix of the distance between segmented volumes i and j as the distance between their centre of mass
    % %the factor 1000 converts the distance from micronmeter to nanometer
    % for i=1:siz_point(1)
    % for j=1:siz_point(1)
    % 
    % distances(i,j)= sqrt( (center_of_mass(i,1)-center_of_mass(j,1))^2 + (center_of_mass(i,2)-center_of_mass(j,2))^2 + (center_of_mass(i,3)-center_of_mass(j,3))^2 )*1000;
    % end
    % end
    % 
    % %Here is determined if two volumes are clustering together i.e. if their distance <= the sum of their radius multiplied for clust_par defined above, or not
    % for i=1:siz_point(1)
    % for j=1:siz_point(1)
    % 
    % if i==j
    %     cluster(i,j)=0;
    % else
    %     if  distances(i,j)<=clast_par*(radius_nn(i)+radius_nn(j))
    %     cluster(i,j)=1;
    %     end
    % end
    % 
    % end
    % end
    % cluster3=cluster; %%%here the cluster matrix is doubled for keep avaliable its information at the end of the script
    % Read the data from Sheet 5 (Center of Homogeneous Mass)
    center_of_mass = readtable(filename, 'Sheet', 5);
    % Display the modified variable names (these are the names used by MATLAB)
    disp('Modified Variable Names:');
    disp(center_of_mass.Properties.VariableNames);
    % Display the original column headers (as stored in the Excel file)
    disp('Original Column Names:');
    disp(center_of_mass.Properties.VariableDescriptions);
    % Define the starting row (based on the repeating headers)
    %start_row = 3; % Adjust this according to your actual data structure
    % Read the table, starting from the correct row
    %center_of_mass = readtable(filename, 'Sheet', 5, 'Range', sprintf('A%d', start_row));
    % Optional: Remove rows that are still unwanted (if any)
    % For example, if there are rows with unwanted headers in the middle of the data:
    %center_of_mass = center_of_mass(~ismember(center_of_mass{:, 1}, {'CenterofHomogeneousMass'}), :);
    %center_of_mass = center_of_mass(~ismember(center_of_mass{:, 1}, {'CenterofHomogeneousMassX'}), :);
    % Use the modified column names to extract the data
    x_coords = center_of_mass{:, 'CenterOfHomogeneousMass'};  % Extract X coordinates
    y_coords = center_of_mass{:, 'Var2'};  % Adjust as needed for Y coordinates
    z_coords = center_of_mass{:, 'Var3'};  % Adjust as needed for Z coordinates
    % Determine the number of detected spots (volumes)
    siz_point = size(center_of_mass);  % Size of the table, number of rows (spots)
    total_resolved_spots = siz_point(1);  % Total number of volumes quantified
    disp(['Total number of volumes quantified: ', num2str(total_resolved_spots)])
    % Create a matrix to store the distances between detected volumes
    disp('Initializing distance matrix...');
    distances = zeros(total_resolved_spots, total_resolved_spots);  % Initialize distance matrix
    % Create a matrix to store cluster information (1 for close volumes, 0 for far volumes)
    disp('Initializing cluster matrix...');
    cluster = zeros(total_resolved_spots, total_resolved_spots);  % Initialize cluster matrix
    % Loop to calculate the distance matrix between detected volumes
    % The factor 1000 converts the distance from micrometers to nanometers
    disp('Calculating distance matrix...');
    for i = 1:total_resolved_spots
        for j = 1:total_resolved_spots
            distances(i,j) = sqrt( (x_coords(i) - x_coords(j))^2 + (y_coords(i) - y_coords(j))^2 + (z_coords(i) - z_coords(j))^2 ) * 1000;  % Calculate distance in nm
        end
    end
    % Loop to calculate whether two volumes belong to the same cluster
    % Check if distance between volumes is less than or equal to the sum of their radii
    disp('Determining cluster assignments...');
    % ERROR_NAME_SHOULD_BE_clast_par.    clust_par = 1.5;  % Define your clustering parameter (this is an example, adjust as necessary)
    for i = 1:total_resolved_spots
        for j = 1:total_resolved_spots
            if i == j
                cluster(i,j) = 0;  % No clustering with itself
            else
                if distances(i,j) <= clast_par * (radius_nn(i) + radius_nn(j))  % Check if the distance is less than the sum of radii
                    cluster(i,j) = 1;  % Assign to the same cluster
                end
            end
        end
    end
    cluster3 = cluster;  % Store the cluster matrix for further use
    %disp('Cluster Matrix (1 indicates volumes in the same cluster, 0 otherwise):');
    %disp(cluster3);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % determination of the number of clusters the cluster size i.e. how many
    % volumes are contained into one cluster and other cluster information
    numclast=0;
    numclast2(1:siz_point(1),1)=1;
    numclast3(1:siz_point(1),1:siz_point(1))=0;
    num=0; % initial number of clusters
    cluster_record(1:siz_point(1),1:siz_point(1))=0;% matrix where the label of the cluster members should be stored i.e the label
    %of the volumes of cluster members. Please note that the column number of cluster_record matrix is also the label of the first volume find for a given
    %cluster if this column is not made only by zero values
    % searching and isolating clusters
    for i=1:siz_point(1)
        count=0;  % how many volumes are clustering with an individual volume
        count3=0; % how many volumes are belonging to the same cluster
        if sum( cluster(i,:))>=1
        num=num+1;
        %     ind= sum( cluster(i,:));
            for j=1:siz_point(1)
                if  cluster(i,j)==1
                    count=count+1;
                    count3=count3+1;
                    cluster(i,j)=0;% every time if a point on the cluster matrix is found to belong to a specific cluster it is set to
                    %zero in order to not choose this point twice. At the end of the cluster search the cluster matrix should be a full zeros matrix
                    cluster(j,i)=0;% everytime if a point on the cluster matrix is found to belong to a specific cluster is set to
                    %zero in order to the choose this point twice at the end the cluster matrix should be a full zeros matrix
                    vect(count)=j;% here are listed the label of the volumes clastering with the volume under analysis
                    cluster_record(count3,i)=j;
                end
            end
        end
        while count > 0;%%%% in this recursive loop are searched the volumes that were clustering with the first one and in case
            %they are found clustering with other volumes these are stored and this cycle is repeated until when no more clustering volume are found
           count2=count;
           count=0;
           vect2=vect;% saving the information contained in vect (is a variable size vector so overwriting is not to suggest)
           clear vect % delete this variable in a way it can be defined newly in the following part of the loop
           for ss=1:count2
                for j=1:siz_point(1)
                if  cluster(vect2(ss),j)==1
                 count=count+1;
                 cluster(vect2(ss),j)=0;
                 cluster(j,vect2(ss))=0;
                vect(count)=j;
                ll=0; %%ll=0 a volume is for the first time found to be member of a given cluster; ll=1 the volume was already found to be part of the cluster
                % and it will be not counted twice
                     for ms=1:count3
                         if cluster_record(ms,i)==j
                             ll=1;
                         end
                     end
                  if ll==0
                count3=count3+1;
                cluster_record(count3,i)=j;
                  end
                  ll=0;
                end
                end
            end
            clear vect2 %  this variable vector will be redefined at the beginning of the loop
            count;
        end
        if num > 0 % if a cluster is found
            if count3>0
                numb(num)=count3+1;% 1 comes from the fact that a cluster is always found starting from a volume that is not included
                %in the counting routine
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dysplaying the information
    sit=size(numb);
    number_of_cluster=sit(2)% the number of columns of the variable numb gives the total number of clusters
    total_number_of_spot_in_cluster= sum(numb)% the total number of volumes belonging to the one of the clusters (it must be inferior to the total number of volumes quantified!!!)
    average_number_of_spot_in_cluster= mean(numb)%  mean number of volumes belonging to a cluster
    dispersion_number_of_spot_in_cluster= std(numb)% variance of the number of volumes belonging to a cluster
    frequent_number_of_spot_in_cluster= median(numb)% median of the number of volumes belonging to a cluster
    indd = 0;
    % Preallocate dist and disti arrays with a sensible size limit
    max_cluster_size = siz_point(1);  % Max number of volumes per cluster, assuming worst case
    % Preallocate clust_info for efficiency
    clust_info = zeros(siz_point(1), 8); % Preallocate space for cluster info
    for l = 1:siz_point(1)
        % Display progress
        fprintf('Processing cluster %d of %d\n', l, siz_point(1));
        % Initialize variables
        cluster_member = [];
        if sum(cluster_record(:, l)) >= 1
            indd = indd + 1;
            cluster_member = [l; cluster_record(cluster_record(:, l) >= 1, l)]; % All labels in the cluster
        end
        if ~isempty(cluster_member)
            got = numel(cluster_member); % Number of volumes in the cluster
            rel = 0; % To keep track of valid distances
            xx = 0; yy = 0; zz = 0;
            % Sum x, y, z positions (in nm, factor 1000)
            cluster_coords = [center_of_mass{cluster_member, 'CenterOfHomogeneousMass'}, ...
                              center_of_mass{cluster_member, 'Var2'}, ...
                              center_of_mass{cluster_member, 'Var3'}] * 1000; % Coordinates in nm
            xx = sum(cluster_coords(:, 1));
            yy = sum(cluster_coords(:, 2));
            zz = sum(cluster_coords(:, 3));
            % Vectorized distance calculation between all points in the cluster
            % We calculate the pairwise distance matrix for all points in the cluster
            dist_matrix = pdist2(cluster_coords, cluster_coords); % pdist2 is highly optimized
            % Extract non-zero distances and store them in disti
            disti = dist_matrix(dist_matrix > 0); % Convert to vector and remove zeros
            % Collect statistical information
            clust_info(indd, 1) = numb(indd); % Number of volumes in the cluster
            clust_info(indd, 2) = 0.5 * mean(disti); % Mean distance between volumes
            clust_info(indd, 3) = std(disti / 2); % Estimate of variance
            clust_info(indd, 4) = max(disti); % Max distance
            clust_info(indd, 5) = min(disti); % Min distance
            clust_info(indd, 6) = xx / got; % X coordinate of the cluster center of mass
            clust_info(indd, 7) = yy / got; % Y coordinate of the cluster center of mass
            clust_info(indd, 8) = zz / got; % Z coordinate of the cluster center of mass
            % Clear temporary variables for next iteration
            clear dist_matrix disti cluster_coords
        end
    end
    % --------- Preparing the final results and save them in the Excel files ---------
    %%% the labels of cluster members are saved in a two column matrix the
    %%% second column contains the label of the first member of the cluster and
    %%% it repeats its value until when the same cluster is under consideration.
    %%% The first column displays the label of the other volumes belonging to
    %%% the clusters
    a=0;
    for i=1:siz_point(1)
        if sum(cluster_record(:,i))>=1
            for j=1:siz_point(1)
                if cluster_record(j,i)>=1
                    a=a+1;
                    position(a,1)=cluster_record(j,i);
                    position(a,2)=i;
                end
            end
        end
    end
    position2= position-1; % since the labelling in Imaris starts from number 0 and not from one in order to find the same cluster in Imaris to the
    % to the label values must be substructed 1
    %
    %%%%%%%% printing in the output file the data  quantified till now
    %%%%%%
    names={'number of particles', 'mean distance', 'variance of distance', 'max distance', 'min distance', 'center_x', 'center_y', 'center_z' };
    names2={'Mean_radius_nm', 'Variance_radius_n ', 'Total_resolved_spots' };
    values=[Mean_radius_nm , Variance_radius_n  total_resolved_spots];
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],names2,1,'A1')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],values,1,'A2')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],clust_info,2,'A2')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],names,2,'A1')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],position2,3,'A1')
    % Write 'names2' and 'values' to the first sheet
    %DEBUG: writetable(table(values), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 1, 'WriteRowNames', false, 'WriteVariableNames', false);
    %DEBUG: writetable(table(names2), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 1, 'WriteRowNames', false);
    % Create a table where the first row is the header and the second is the data
    output_table = [names2; num2cell(values)];
    disp(output_table);
    % Write the table to the Excel file, starting at cell A1
    writetable(cell2table(output_table), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 1, 'WriteRowNames', false, 'WriteVariableNames', false);
    % % Saving as an older .xls format (Excel 97-2003 format)
    % writetable(cell2table(output_table), [num2str(clast_par) '_output_' filename '.xls'], 'Sheet', 1, 'WriteRowNames', false, 'WriteVariableNames', false);
    % Save as CSV file
    writetable(cell2table(output_table), [num2str(clast_par) '_output_' filename '_Sheet1.csv'], 'WriteRowNames', false, 'WriteVariableNames', false);
    % Convert 'clust_info' into a table if it's not one already and write to the second sheet
    clust_info_table = array2table(clust_info, 'VariableNames', names);  % Convert 'clust_info' to table with proper column names
    writetable(clust_info_table, [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 2, 'WriteRowNames', false);
    % Assuming position2 is a matrix that needs to be written to the third sheet
    writetable(array2table(position2), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 3, 'WriteRowNames', false, 'WriteVariableNames', false);
    %%%% quantification of the distance between different clusters
    %%%%% the distance between clusters is defined to be the distance between their centre of
    %%%%% mass
    dist_clust(1:sit(2),1:sit(2))=0;
    cc1=0;
    for i=1:sit(2)
        for j=1:sit(2)
            if i==j
                dist_clust(i,j)=100000;% the fact that the distance of one cluster respect itself is zero would make more complicate
                %to determine the minimum distance between a cluster and it neighbours
            else
                cc1=cc1+1;
                dist_clust(i,j)=sqrt( (clust_info(i,6)-clust_info(j,6))^2 + (clust_info(i,7)-clust_info(j,7))^2 + (clust_info(i,8)-clust_info(j,8))^2 );
                dist_clust2(cc1)=dist_clust(i,j);
            end
        end
    end
    mindist_cluster=min(dist_clust);% the minimum distance between one cluster and the other ones
    mean_dist_cluster=mean(mindist_cluster);% the mean distance between closer clusters
    %%% exporting the latest information
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'], dist_clust,4,'A1')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],mean_dist_cluster,5,'A1')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],mindist_cluster',5,'B1')
    % names3={'number_of_cluster', 'total_number_of_spot_in_cluster ', 'average_number_of_spot_in_cluster', 'dispersion_number_of_spot_in_cluster','frequent_number_of_spot_in_cluster' }
    % values3=[number_of_cluster total_number_of_spot_in_cluster average_number_of_spot_in_cluster dispersion_number_of_spot_in_cluster frequent_number_of_spot_in_cluster];
    % values3=values3';
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],names3,6,'A1')
    % xlswrite([num2str(clast_par) '_output_' filename '.xlsx'],values3',6,'A2')
    % Exporting the latest information
    writetable(table(dist_clust), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 4, 'WriteRowNames', true);
    writetable(table(mean_dist_cluster), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 5, 'WriteRowNames', true);
    writetable(table(mindist_cluster'), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 5, 'WriteRowNames', true, 'WriteVariableNames', false);
    % Save other information about clusters
    names3 = {'number_of_cluster', 'total_number_of_spot_in_cluster', 'average_number_of_spot_in_cluster', ...
        'dispersion_number_of_spot_in_cluster', 'frequent_number_of_spot_in_cluster'};
    values3 = [number_of_cluster, total_number_of_spot_in_cluster, average_number_of_spot_in_cluster, ...
        dispersion_number_of_spot_in_cluster, frequent_number_of_spot_in_cluster];
    %DEBUG: writetable(table(values3), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 6, 'WriteRowNames', true, 'WriteVariableNames', false);
    %DEBUG: writetable(table(names3), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 6, 'WriteRowNames', true);
    % Create a table where the first row is the header and the second row is the data
    output_table3 = [names3; num2cell(values3)];
    disp(output_table3);
    % Write the combined table to the Excel file, starting from cell A1 of Sheet 6
    writetable(cell2table(output_table3), [num2str(clast_par) '_output_' filename '.xlsx'], 'Sheet', 6, 'WriteRowNames', false, 'WriteVariableNames', false);
    % % Saving as an older .xls format (Excel 97-2003 format)
    % writetable(cell2table(output_table3), [num2str(clast_par) '_output_' filename '.xls'], 'Sheet', 6, 'WriteRowNames', false, 'WriteVariableNames', false);
    % Save as CSV file
    writetable(cell2table(output_table3), [num2str(clast_par) '_output_' filename '_Sheet6.csv'], 'WriteRowNames', false, 'WriteVariableNames', false);