%% Antibody Nuclear and Cytosolic Intensity Quantification

% Authors: Dr. Joshua Morgan and Isabella Bagdasarian, TIME lab, UCR
% use this script to get "quick" IF intensity quantifications of nuclear and
% cytosolic antibodies. It does not necessarily run fast (that depends on file size),
% by quick I mean it does not require much thought or organization to run this.

%%% Outputs a folder containing traced images, a folder of antibody images
%%% and a .mat file containing the name of each series, nuclear antibody intensity,
%%% cytosolic antibody intensity, and the ratio of nuclear to cytosolic
%%% antibody intensity.

%%% March 17th, 2020: updated to process multiple lif files at once

clearvars
clc
close all

tEntire = tic;

%image output path
outdir1 = 'C:\Users\jorda\Desktop\TIME Lab\ASC Mitochondrial Morphology Experiment\Intensity Quantification\';


%change these things
IM_FILE_all = ["20231026_ASC_GFP_MitoLocal_MAHS_T4_DAPI_GFP_Phalloidin554_MitoView640_TS6x6_20x_JLR.lif"];

antibody = 'MitoView'; %label the antibody you want to quantify

%change this if segmentation is not working

%%%%greedy nuclear segmentation parameters
gNthresh = 10;
gNAreaOpen = 100; %was 350 was 100

%%%%strict nuclear segmentation parameters
sNAreaOpen = 50; %strict nuclear area open %was 150 was 50

%%%%phalloidin segmentation parameters
gPthresh = 1000; %changed from 6000 to 1400 on March17th2020 %changed to 1000 on Jan30th2024
gPAreaOpen = 600; 

%Set to true to generate images
wantImages = true;


%length(IM_FILE_all)
for file = 1:length(IM_FILE_all)
    
    IM_FILE = (sprintf('%s',IM_FILE_all(file,1)));
    
    %open and setup output folder
    outdir2 = sprintf("%s_%s",IM_FILE(1:end-4), 'data output'); %create a folder within analysis directory
    
    if ~isfolder(outdir2)
        mkdir(outdir2);
    end
    
    %open and setup output for traced nuclear and cytosolic images
    if ~isfolder(sprintf('%s\\%s\\traced image output %s',outdir1, outdir2,IM_FILE(1:end-4)))
        mkdir(sprintf('%s\\%s\\traced image output %s',outdir1, outdir2,IM_FILE(1:end-4)));
    end
    
    %open and setup output for antibody images
    if ~isfolder(sprintf('%s\\%s\\%s image output %s',outdir1, outdir2,antibody,IM_FILE(1:end-4)))
        mkdir(sprintf('%s\\%s\\%s image output %s',outdir1, outdir2,antibody,IM_FILE(1:end-4)));
    end
    
    tic
    rdr = bfGetReader(IM_FILE); %the reader variable is a pointer to the file
    omeMeta = rdr.getMetadataStore(); %pull OME metadata
    nS = rdr.getSeriesCount(); %number of series
    toc
    
    %preallocate
    Sname = cell(nS,1);
    intDataNuc = cell(nS,1); %for nuclear intensity data
    intDataCyt = cell(nS,1); %for cytosolic intensity data
    intDataRat = cell(nS,1); %for nuclear/cytosolic ratio intensity data
    
    tic;A = bfopen(IM_FILE);toc; %reads file
    fprintf('finished opening file');
    
    for g = 1:length(A)
        
        %collect information about the series
        %since JAVA objects (like the OME metadata) start there indexing at 0, we
        %set the reader to be 0.
        rdr.setSeries(g-1);
        nC = omeMeta.getChannelCount(g-1); %number of channels
        nP = omeMeta.getPlaneCount(g-1); %number of planes
        SZ = omeMeta.getPixelsSizeZ(g-1).getValue(); % number of Z slices
        SX = omeMeta.getPixelsSizeX(g-1).getValue(); % image width, pixels
        SY = omeMeta.getPixelsSizeY(g-1).getValue(); % image height, pixels
        Sname{g} = char(omeMeta.getImageName(g-1)); %image name
        iP = nP/nC;
        series_1 = g; %current series
        
        %preallocate arrays
        D1 = zeros(SY,SX,'uint16');
        Y1 = zeros(SY,SX,'uint16');
        P1 = zeros(SY,SX,'uint16');
        
        %create a max projection by sequentially comparing the current image to
        %the next slice and overwriting pixels based on max
        
        D1 = cat(3,A{series_1,1}{1:iP,1}); %ch1 %nuclei channel
        Y1 = cat(3,A{series_1,1}{iP*3+1:iP*4,1}); %ch2 %antibody channel (MitoView)
        P1 = cat(3,A{series_1,1}{iP*2+1:iP*3,1}); %ch3 %f-actin/cell body marker channel
        
        %use focus measure to find a maximum projection
        [D1] = laplacianfocusMAX(D1);
        [Y1] = laplacianfocusMAX(Y1);
        [P1] = laplacianfocusMAX(P1);
        
        %% DRAQ7 segmentation
        
        [Do] = NucCleanUp(D1);
        
        %use a nuclei by nuclei threshold approach
        %perform initial greedy segmentation
        [Db2] = GreedyNucSeg(Do,  gNthresh, gNAreaOpen);
        
        %perform second stricter segmentation
        [Db4, ~, Df] = StrictNucSeg(Db2, Do, sNAreaOpen);
        
        [Dws] = NucWatershed_2D(Db4);
        
        %Peform one last step to create a shrunk version for the cytosolic step
        %below as well as a thicken version with the same connectivity
        Db4e = bwmorph(Dws,'fill'); %#ok<*NASGU>
        Db4e = bwmorph(Dws,'shrink',2);
        Db4et = bwmorph(Dws,'thicken',2); %thicken without combining discconnected objects
        
        %Reject phalloidan artifacts
        [Ld] = CancelCytArt(D1, Db4et);
        
        %% Phalloidan segmentation
        [Po, Pth] = CytCleanUp(P1);
        
        
        [T, mask] = CytGradSeg(Pth, Po, Db4e, gPAreaOpen, gPthresh);
        L = watershed(T);
        L(~mask) = 0;
        Lp = L;
        Lp_org = Lp;
        
        %% Fix and correlate Nuc and Cyt labels
        
        %delete & renumber cytosol that have no nuclei (likely artifact)
        TEST1 = uint8(Ld>0)+uint8(Lp>0); %create image that is 1 for cytosol and 2 for cytosol+nuclei
        TEST2 = hysthresh(TEST1,0.5,1.5); %use hysteresis thresholding to eliminate cytosol without nuclei
        Lp(TEST2==0) = 0; %delete
        CC = bwconncomp(Lp>0); %recreate connected components
        Lp = labelmatrix(CC); %renumber nuclear label matrix
        
        %delete & renumber nuclei that have no cytosol (likely artifact)
        Ld(Lp==0) = 0; %delete
        Ld = bwareaopen(Ld,100); %delete "fragments" caused by irregular nuc/cyto overlays
        CC = bwconncomp(Ld>0); %recreate connected components
        Ld = labelmatrix(CC); %renumber nuclear label matrix
        
        %Match cytosol to nuclei
        %first, erode the Ld label matrix (from Db4et) down to Db4e size
        Lde = Ld;
        Lde(~Db4e) = 0;
        CtN = regionprops(Lp,Lde,'MaxIntensity');
        CtN = [CtN.MaxIntensity];
        Lp_new = zeros(size(Lp),class(Lp)); %create a new (blank) cytosolic label matrix
        for CtNi = 1:length(CtN) %loop over each entry in CtN and change the Lp patch to the new number given in CtN
            Lp_new(Lp == CtNi) = CtN(CtNi);
        end
        TEST = Lp_new; TEST(Ld>0) = Ld(Ld>0);
        
        %you can use the below line to test. If you see nuclei, it isn't well
        %matched up
        %figure; image(TEST,'cdatam','scal');colormap('colorcube')
        
        %delete from masks
        Db4et_org = Db4et;
        Db4et(Ld==0) = 0;
        
        %% Create Boundaries
        
        %Identify Draq7 Boundaries
        B = bwboundaries(Db4et, 'noholes');
        
        %identify phalloidin boundaries
        B2 = cell(max(Lp(:)),1);
        for l = 1:max(Lp(:))
            B2{l} = bwboundaries(Lp==l, 'noholes');
        end
        
        %Identify Draq7 Boundaries - orginal
        Bo = bwboundaries(Db4et_org, 'noholes');
        
        %identify phalloidin boundaries
        B2o = cell(max(Lp_org(:)),1);
        for l = 1:max(Lp_org(:))
            B2o{l} = bwboundaries(Lp_org==l, 'noholes');
        end
        
        %Calculate antibody intensity in the nuclei
        SYn = regionprops(Ld, Y1, 'MeanIntensity','Centroid');
        centroid = cat(1,SYn.Centroid);
        Y_int = [SYn.MeanIntensity]';
        intDataNuc{g} = Y_int;
        
        %Calculate antibody cytosolic intensity
        Lc = Lp_new;
        Lc(Db4et) = 0; %remove nuclei from the cytosol, so we only measure cyt.
        SYc = regionprops(Lc, Y1, 'MeanIntensity');
        C_int = [SYc.MeanIntensity]';
        intDataCyt{g} = C_int;
        
        %ratio of nuclear to cytosolic intensity
        if size(Y_int) == size (C_int)
        intDataRat{g} = Y_int ./ C_int;
        else
            continue
        end
        toc
        
        if wantImages
            IM = cat(3,imadjust(P1),zeros(size(Y1)), imadjust(D1));
            image(IM,'cdatamapping','scaled')
            hold on
            
            %Plot trace of phalloidin
            for s = 1:size(B2,1)
                UN = B2{s,1};
                for f2 = 1:length(UN)
                    [boundary2] = UN{f2};
                    plot(boundary2(:,2), boundary2(:,1), 'c', 'LineWidth', 1)
                end
            end
            
            %Plot trace of nuclei
            for f = 1:length(B)
                boundary = B{f};
                plot(boundary(:,2), boundary(:,1), 'm', 'LineWidth', 1)
            end
            
            %uncomment the below section to see the original tracings (before
            %correlation correction)
            %         %Plot trace of phalloidin-original
            %         for s = 1:size(B2o,1)
            %             UN = B2o{s,1};
            %             for f2 = 1:length(UN)
            %                 [boundary2] = UN{f2};
            %                 plot(boundary2(:,2), boundary2(:,1), 'c:', 'LineWidth', 2)
            %             end
            %         end
            %
            %         %Plot trace of nuclei-original
            %         for f = 1:length(Bo)
            %             boundary = Bo{f};
            %             plot(boundary(:,2), boundary(:,1), 'm:', 'LineWidth', 2)
            %         end
            hold off
            drawnow
            
            print(sprintf('%s\\%s\\traced image output %s\\series%d' , outdir1, outdir2,IM_FILE(1:end-4),g),'-djpeg','-r300');
            
            % Save antibody channel
            figure(2)
            IM = cat(3,zeros(size(D1)), imadjust(Y1), zeros(size(D1)));
            image(IM,'cdatamapping','scaled'); movegui('northeast')
            hold off
            drawnow
            print(sprintf('%s\\%s\\%s image output %s\\series%d' , outdir1, outdir2,antibody, IM_FILE(1:end-4),g),'-djpeg','-r300');
            
        end
        close all
        
    end
    save(sprintf('%s\\%s\\%s_data.mat',outdir1,outdir2,IM_FILE(1:end-4)),'Sname','intDataNuc','intDataCyt','intDataRat')
    
   % clearvars -except tEntire IM_FILE_all outdir1 outdir2 file wantImages antibody gNthresh gNAreaOpen ...
  %     sNAreaOpen gPthresh gPAreaOpen
end

tCompEnt = toc(tEntire);
fprintf(1,'Analysis took %d hours to complete.\n',tCompEnt/3600);

toc;