/* * Erosion Territories * * Expected Input: * The expected input for this script are 2D non-timelapse images in tiff/tif/ND2 format that are either 3 or 4 channel * One channel is expect to be a DAPI or other nuclear stain * At start of macro, user selects folder with images in (nested folders dealth with fine), sets channel order and chooses desired number of bins * * Description of script workflow: * The DAPI channel is used to segment the nuclei (using gaussian blur and huang threshold and the 'watershed irregular features' from BioVoxxel - see installation notes) * Each nucleus area is duplicated from original image to isolate nuclei for further processing * For each nuclei image, the area is measured and the expected area of the bins is calculated based on the number requsted * Using the expected areas of these bins and iteratively thresholding on a distance map of the nucleus ROI bins are produced * Mean intensity of these ROIs is measured in all the channels * * Output: * Each nucleus found in each image will be saved as an RGB image with the outline of created bins overlaid on the image * A .csv file will be saved with the areas of the nuclei, the bins and the mean intensity of each bin for each channel in the image * Bin/ territory naming in the results go from outside in so bin 1 is the outermost one. * * Installation Notes: * Download FIJI here - https://imagej.net/software/fiji/downloads * Details of the Biovoxxel Update site - https://imagej.net/plugins/biovoxxel-toolbox * How to add an update site - https://imagej.net/update-sites/ * How to use an ImageJ macro from Github - https://github.com/IGC-Advanced-Imaging-Resource/ImageJ_Macros * This script is hosted here with an example image - https://github.com/IGC-Advanced-Imaging-Resource/Purshouse2022_paper * * Written by Laura Murphy (laura.murphy@ed.ac.uk) * IGC Advanced Imaging Resource (https://www.ed.ac.uk/institute-genetics-cancer/facilities/advanced-imaging-resource) * First written: November 2019 Last updated January 2022 */ //--------------------------------//------------------------------------------------------------------------------------------------------- //-- Part 0: Preparation steps: get directories from user, set up arrays to store results that are being added to and some other setup //--------------------------------//------------------------------------------------------------------------------------------------------- // Get user to select folder with images for input and the output folder where they want results saved input_folder = getDirectory("Parent directory of your image folders"); output_folder = getDirectory("Choose the folder where you want to save your results"); // Set up empty arrays that will be filled during script and displayed and saved at the end Image_Name = newArray(); Nucleus_No = newArray(); Nucleus_Area = newArray(); Bin_No = newArray() Bin_Area = newArray() DAPI_MeanIntensity = newArray(); FITC_MeanIntensity = newArray(); TxRd_MeanIntensity = newArray(); Cy5_MeanIntensity = newArray(); // Magic numbers for nuclei segmentation filter_sigma = 2; // Sigma for Gaussian blur erosion_cycle = 1; // Erosion cycle for watershed watershed_convexity = 0.99; // Convexity threshold for watershed min_particle_size = 50; // Minimum particle size for Analyze Particles max_particle_size = 700; // Maxmimum particle size for Analyze Particles min_particle_circularity = 0.7; // Minimum circularity for Analyze Particles max_particle_circularity = 1.00; // Maxmimum circularity for Analyze Particles // Array of channel colours colour_array = newArray("Blue", "Green", "Red", "Cyan"); roiManager("Reset"); // Get list of images in nested folders with certain extensions - function at bottom of macro directory_list = newArray(); directory_list = get_file_tree(input_folder, directory_list); //--------------------------------//-------------------------------------------------------------------------------------------- //-- Part 1: Image loop starts, each image goes through nuclear segmentation and each nuclei is duplicated from original image //--------------------------------//-------------------------------------------------------------------------------------------- // Start looping through images in input folders for (img = 0; img < directory_list.length; img++) { path = directory_list[img]; run("Bio-Formats Importer", "open=[" + path + "] autoscale color_mode=Grayscale rois_import=[ROI manager] view=Hyperstack stack_order=XYCZT"); getDimensions(width, height, channels, slices, frames); original_image = getTitle(); image_name = File.nameWithoutExtension; // When first image is open, get user to input desired number of bins and assign channels in order of the stack if (img == 0) { Dialog.create("User input"); Dialog.addMessage("Enter the number of concentric territories you want the nucleus split into"); Dialog.addNumber("Bins", 4); Dialog.addMessage("Enter the order of your channels below"); Dialog.addNumber("DAPI", 4); Dialog.addNumber("FITC", 1); Dialog.addNumber("TxRd", 2); Dialog.addNumber("Cy5", 3); Dialog.show; //Store user input bins = Dialog.getNumber; channel_order = newArray(channels); for (ch = 0; ch < channels; ch++){ channel_order[ch] = Dialog.getNumber; } } // Get pixel sizing, if image isn't scaled then skip image getVoxelSize(width, height, depth, unit); if (unit == "pixels"){ waitForUser("Image has no scaling and will be skipped"); run("Close"); continue; } // Set channel colours for (lut = 0; lut < channel_order.length; lut++){ selectWindow(original_image); Stack.setChannel(channel_order[lut]); run(colour_array[lut]); } // Duplicate image to preserve for display use selectWindow(original_image); run("Duplicate...", "title=Display duplicate"); // Segment nuclei in image selectWindow(original_image); Stack.setChannel(channel_order[0]); run("Duplicate...", "title=Nuclear_Mask"); run("Gaussian Blur...", "sigma=" + filter_sigma); setAutoThreshold("Huang dark"); run("Convert to Mask"); run("Watershed Irregular Features", "erosion=" + erosion_cycle + " convexity_threshold=" + watershed_convexity + " separator_size=0-Infinity"); run("Analyze Particles...", "size=" + min_particle_size + "-" + max_particle_size + " circularity=" + min_particle_circularity + "-" + max_particle_circularity + " exclude include add"); // Store number of nuclei number_of_nuclei = roiManager("Count"); // Loop through nuclei, duplicate image of each one in image, naming them so they can be processed further for (nuclei = 0; nuclei < number_of_nuclei; nuclei++){ selectWindow("Nuclear_Mask"); roiManager("Select", nuclei); title = "Nuclei_" + nuclei + 1; run("Duplicate...", "title=" + title); selectWindow(original_image); roiManager("Select", nuclei); run("Duplicate...", "duplicate"); rename("Measure" + title); selectWindow("Display"); roiManager("Select", nuclei); run("Duplicate...", "duplicate"); rename("Display" + title); } //--------------------------------//-------------------------------------------------------------------------------------------- //-- Part 2: Looping through isolated nuclei image starts, creation of bins //--------------------------------//-------------------------------------------------------------------------------------------- // Start of loop for going through the duplicated images of the solo nuclei for (image_of_nucleus = 0; image_of_nucleus < number_of_nuclei; image_of_nucleus++){ title = "Nuclei_" + image_of_nucleus + 1; roiManager("Reset"); selectWindow(title); // Create distance map to use for erosion run("Distance Map"); run("Clear Outside"); // Get area of whole nuclei and maximum value of distance map (further distance in nuclei mask from edge) getStatistics(whole_area, mean, min, whole_max, std, histogram); nuclear_area = whole_area; run("Select None"); // Work out the desired area each bin should be by dividing total nuclear area by number of bins area_fraction = whole_area/bins; prebins = newArray(bins); // Work out the expected areas of the shells (i.e. the smaller circles that share centroid with nucleus) for (areas_of_prebins = 0; areas_of_prebins < prebins.length; areas_of_prebins++){ addtoArray = whole_area; prebins[areas_of_prebins] = addtoArray; whole_area = whole_area-area_fraction; } // For each of the chosen bin number, iteratively reduce threshold of distance map until area of thresholded area is lower/equal to expected // Once that is reached, add to ROI Manager. for (bin = 0; bin < prebins.length; bin++) { current_area_aim = prebins[bin]; for (threshold = 1; threshold < whole_max; threshold++) { getStatistics(area, mean, min, max, std, histogram); run("Select None"); setThreshold(threshold, whole_max); run("Create Selection"); getStatistics(area, mean, min, max, std, histogram); if (area <= current_area_aim){ roiManager("Add"); break; continue; } } } // Use prebins to create bins by using XOR function between the pre-bins in the ROI Manager number_of_rois = roiManager("Count"); to_be_deleted = newArray(number_of_rois); // Add current ROIs to array to delete later for (bins = 0; bins < number_of_rois; bins++){ to_be_deleted[bins] = bins; } // Loop through rois, selecting pairs to get XOR area between them, add to ROI and name them. for (roi = 0; roi < number_of_rois-1; roi++){ roiManager("select", newArray(roi, roi+1)); roiManager("XOR"); Roi.setName("BIN_" + roi+1); roiManager("Add"); roiManager("select", newArray(roi, roi+1)); } roiManager("select", number_of_rois-1); Roi.setName("BIN_" + number_of_rois); roiManager("Add"); // Delete the shells from earlier roiManager("select", to_be_deleted); roiManager("Delete"); resetThreshold(); //--------------------------------//-------------------------------------------------------------------------------------------- //-- Part 3: Time to save results, for each bin, the image name, nucleus and nucleus area is saved to keep arrays same length // Then for each of the bins (also called territories), the area is measured and then the intensity of each of the channels // All things are saved to arrays which are displayed once complete. On the last image, the arrays are saved to the output folder //--------------------------------//-------------------------------------------------------------------------------------------- // Loop through final shapes, adding nuclei results to arrays to keep arrays even for (territories = 0; territories < roiManager("Count"); territories++){ Image_Name = Array.concat(Image_Name, image_name); Nucleus_No = Array.concat(Nucleus_No, image_of_nucleus+1); Nucleus_Area = Array.concat(Nucleus_Area, nuclear_area); Bin_No = Array.concat(Bin_No, territories+1); selectWindow("Measure" + title); roiManager("Select", territories); getStatistics(bin_area, mean, min, max, std, histogram); Bin_Area = Array.concat(Bin_Area, bin_area); Stack.setChannel(channel_order[0]); roiManager("Select", territories); getStatistics(bin_area, mean_DAPI, min, max, std, histogram); DAPI_MeanIntensity = Array.concat(DAPI_MeanIntensity, mean_DAPI); Stack.setChannel(channel_order[1]); roiManager("Select", territories); getStatistics(bin_area, mean_FITC, min, max, std, histogram); FITC_MeanIntensity = Array.concat(FITC_MeanIntensity, mean_FITC); Stack.setChannel(channel_order[2]); roiManager("Select", territories); getStatistics(binArea, mean_TXRD, min, max, std, histogram); TxRd_MeanIntensity = Array.concat(TxRd_MeanIntensity, mean_TXRD); // If the images are 4 channel, then save the Cy5 channel if (channel_order.length > 3){ Stack.setChannel(channel_order[3]); roiManager("Select", territories); getStatistics(binArea, mean_Cy5, min, max, std, histogram); Cy5_MeanIntensity = Array.concat(Cy5_MeanIntensity, mean_Cy5); } } // Take original image, convert to RGB and save with bin boundaries overlayed on image for checking. selectWindow("Display" + title); run("Make Composite"); run("RGB Color"); roiManager("Show all without labels"); run("Flatten"); saveAs("Tiff", output_folder + image_name + "_" + title); } // Set up for next image coming up by closing open images and reset ROI manager and results table run("Close All"); roiManager("Reset"); run("Clear Results"); // Show updated array of results each time an image is completed. If image was 3 channel, the Cy5 array will be empty. Array.show(Image_Name, Nucleus_No, Nucleus_Area, Bin_No, Bin_Area, DAPI_MeanIntensity, FITC_MeanIntensity, TxRd_MeanIntensity, Cy5_MeanIntensity); } // Save table of arrays to user directed output folder now all images are processed saveAs("Results", output_folder + File.separator + "Results.csv"); run("Close"); //--------------------------------//----------------------------------------------------------------------- //-- Part 4: Let user know the macro is complete //--------------------------------//----------------------------------------------------------------------- Dialog.create("Progress"); Dialog.addMessage("Macro Complete!"); Dialog.show; //--------------------------------//----------------------------------------------------------------------- //-- Functions //--------------------------------//----------------------------------------------------------------------- function get_file_tree(directory, file_tree){ file_list = getFileList(directory); for(f = 0; f < file_list.length; f++){ if (matches(file_list[f], "(?i).*\\.(tif|tiff|nd2)$")) file_tree = Array.concat(file_tree, directory + file_list[f]); if(File.isDirectory(directory + File.separator + file_list[f])) file_tree = get_file_tree(directory + file_list[f], file_tree); } return file_tree;