//******* // Author: Pradeep Rajasekhar // March 2021 // License: BSD3 // // Copyright 2021 Pradeep Rajasekhar, Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia // // Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. // 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. var fs=File.separator; setOption("ExpandableArrays", true); print("\\Clear"); var fiji_dir=getDirectory("imagej"); var gat_dir=fiji_dir+"scripts"+fs+"GAT"+fs+"Tools"+fs+"commands"; //settings for GAT gat_settings_path=gat_dir+fs+"gat_settings.ijm"; if(!File.exists(gat_settings_path)) exit("Cannot find settings file. Check: "+gat_settings_path); run("Results... ", "open="+gat_settings_path); //training_pixel_size=parseFloat(Table.get("Values", 0)); //0.7; neuron_area_limit=parseFloat(Table.get("Values", 1)); //1500 neuron_seg_lower_limit=parseFloat(Table.get("Values", 2)); //90 neuron_lower_limit=parseFloat(Table.get("Values", 3)); //160 backgrnd_radius=parseFloat(Table.get("Values", 4)); run("Close"); //check if required plugins are installed var check_plugin=gat_dir+fs+"check_plugin.ijm"; if(!File.exists(check_plugin)) exit("Cannot find check plugin macro. Returning: "+check_plugin); runMacro(check_plugin); //check if label to roi macro is present var label_to_roi=gat_dir+fs+"Convert_Label_to_ROIs.ijm"; if(!File.exists(label_to_roi)) exit("Cannot find label to roi script. Returning: "+label_to_roi); //check if roi to label macro is present var roi_to_label=gat_dir+fs+"Convert_ROI_to_Labels.ijm"; if(!File.exists(roi_to_label)) exit("Cannot find roi to label script. Returning: "+roi_to_label); //check if ganglia cell count is present var ganglia_cell_count=gat_dir+fs+"Calculate_Neurons_per_Ganglia.ijm"; if(!File.exists(ganglia_cell_count)) exit("Cannot find ganglia cell count script. Returning: "+ganglia_cell_count); //check if ganglia prediction post processing macro present var deepimagej_post_processing=gat_dir+fs+"Ganglia_prediction_post_processing.ijm"; if(!File.exists(deepimagej_post_processing)) exit("Cannot find roi to label script. Returning: "+deepimagej_post_processing); //check if ganglia prediction post processing macro present var segment_ganglia=gat_dir+fs+"Segment_Ganglia.ijm"; if(!File.exists(segment_ganglia)) exit("Cannot find segment ganglia script. Returning: "+segment_ganglia); #@ File (style="open", label="Choose the image to segment.
Enter NA if image is open.") path #@ boolean image_already_open #@ String(value="If image is already open, tick above box.", visibility="MESSAGE") hint1 #@ File (style="open", label="Choose the StarDist model file if segmenting neurons.
Enter NA if empty",value="NA", description="Enter NA if nothing") neuron_model_path #@ String(label="Enter channel number for Hu if you know. Leave as NA if not using.", value="NA") cell_channel #@ String(value="-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------",visibility="MESSAGE") hint_star #@ String(value="
NEURONAL SUBTYPE ANALYSIS
",visibility="MESSAGE") hint_subtype #@ boolean Calculate_Neuron_Subtype #@ String(value="Tick above box if you want to estimate proportion of neuronal subtypes.", visibility="MESSAGE") hint3 #@ File (style="open", label="Choose the StarDist model for subtype segmentation.
Enter NA if empty",value="NA", description="Enter NA if nothing") subtype_model_path cell_type="Neuron"; #@ String(value="If you already know the channel names and numbers, check the box below and enter them.
The channel numbers MUST match the channel name order.
You have the option of entering them later in the analysis",visibility="MESSAGE") hint5 #@ boolean Enter_channel_details_now #@ String(label="Enter channel names followed by a comma (,). Leave as NA if not using.", value="NA") marker_names_manual #@ String(label="Enter channel numbers with separated by a comma (,). Leave as NA if not using.", value="NA") marker_no_manual #@ String(value="-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------",visibility="MESSAGE") hint_star #@ String(value="
DETERMINE GANGLIA OUTLINE
",visibility="MESSAGE") hint_ganglia #@ String(value=" You will get cell counts for each ganglia
If you have a channel for neuron and another channel that labels the ganglia (PGP9.5/GFAP/NOS/Calbindin...)
it should be sufficient to calculate ganglia outline. You can also manually draw the ganglia",visibility="MESSAGE") hint4 #@ String(label=" Enter the channel to use for segmenting ganglia.
Preferably a bright marker that labels most of the ganglia.
Leave as NA if not using. ", value="NA") ganglia_channel #@ boolean Cell_counts_per_ganglia (description="Use a pretrained deepImageJ model to predict ganglia outline") #@ String(choices={"DeepImageJ","Manually draw ganglia"}, style="radioButtonHorizontal") Ganglia_detection #@ String(value="-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------",visibility="MESSAGE") hint_star #@ String(value="
Specify custom scaling (Default is pixel size of 0.7 microns)
",visibility="MESSAGE") hint_scaling #@ String(value="Choose either XY pixel size (microns) or scaling factor (scales images by the specified factor)", visibility="MESSAGE") hint_scaling1 #@ String(choices={"Use pixel size","Use scaling factor"}, style="radioButtonHorizontal",label="Choose scaling option") scaling_option #@ Double (label="Enter scaling value", value=0.568) scale_factor_1 #@ String(value="-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------",visibility="MESSAGE") hint_star #@ String(value="
Finetune cell detection
",visibility="MESSAGE") hint_stardist #@ String(value="Probability determines the number of cells, low values will give more cells.
Reducing overlap threshold will lead to more overlapping cells.
More info about below parameters can be found here: https://www.imagej.net/StarDist/",visibility="MESSAGE", required=false) hint34 #@ String(value="Neuron detection",visibility="MESSAGE") hint_stardist1 #@ Double (label="Probability ", style="slider", min=0, max=1, stepSize=0.05,value=0.55) probability #@ Double (label="Overlap Threhshold", style="slider", min=0, max=1, stepSize=0.05,value=0.5) overlap #@ String(value="Detection of neuronal subtypes",visibility="MESSAGE") hint_stardist2 #@ Double (label="Probability ", style="slider", min=0, max=1, stepSize=0.05,value=0.55) probability_subtype #@ Double (label="Overlap Threhshold", style="slider", min=0, max=1, stepSize=0.05,value=0.5) overlap_subtype //add an option for defining a custom scaling factor marker_subtype=Calculate_Neuron_Subtype; //checking if no of markers and no of channels match if(marker_subtype==1 && Enter_channel_details_now==1) { marker_names_manual=split(marker_names_manual, ","); marker_no_manual=split(marker_no_manual, ","); if(marker_names_manual.length!=marker_no_manual.length) exit("Number of marker names and marker channels do not match"); } if(image_already_open==true) { waitForUser("Select Image. and choose output folder in next prompt"); file_name=getTitle(); //get file name without extension (.lif) dir=getDirectory("Choose Output Folder"); //file_name=File.nameWithoutExtension; } else { if(endsWith(path, ".czi")|| endsWith(path, ".lif")) run("Bio-Formats", "open=["+path+"] color_mode=Composite rois_import=[ROI manager] view=Hyperstack stack_order=XYCZT"); else if (endsWith(path, ".tif")|| endsWith(path, ".tiff")) open(path); else exit("File type not recognised. Tif, Lif and CZI files supported."); dir=File.directory; file_name=File.nameWithoutExtension; //get file name without extension (.lif) } file_name_length=lengthOf(file_name); if(file_name_length>50) file_name=substring(file_name, 0, 39); //Restricting file name length as in Windows long path names can cause errors img_name=getTitle(); Stack.getDimensions(width, height, sizeC, sizeZ, frames); run("Select None"); run("Remove Overlay"); getPixelSize(unit, pixelWidth, pixelHeight); if(unit!="microns") exit("Image not calibrated"); if(scaling_option=="Use pixel size") { //training_pixel_size=scale_factor_1;//0.568; //Images were trained in StarDist using images of this pixel size. Change this for adult human. ~0.9? //Training images were pixelsize of ~0.568, so scaling images based on this scale_factor=pixelWidth/scale_factor_1; if(scale_factor<1.001 && scale_factor>1) scale_factor=1; } else { scale_factor=scale_factor_1; } print(scale_factor); print("Analysing: "+file_name); analysis_dir= dir+"Analysis"+fs; if (!File.exists(analysis_dir)) File.makeDirectory(analysis_dir); print("Files will be saved at: "+analysis_dir); print("Analysing: "+file_name); //Create results directory with file name in "analysis" results_dir=analysis_dir+file_name+fs; //directory to save images if (!File.exists(results_dir)) File.makeDirectory(results_dir); //create directory to save results file //do not include cells greater than 1500 micron in area //neuron_area_limit=1500; //microns neuron_max_pixels=neuron_area_limit/pixelWidth; //convert micron to pixels //do not include cells lower than 90 micron in area for marker //neuron_seg_lower_limit=90;//microns neuron_seg_lower_limit=neuron_seg_lower_limit/pixelWidth; //using limit for marker segmentation //neuron_lower_limit= 160;//microns neuron_min_pixels=neuron_lower_limit/pixelWidth; //convert micron to pixels table_name="Analysis_"+cell_type+"_"+file_name; Table.create(table_name);//Final Results Table row=0; //row counter for the table image_counter=0; if(cell_channel!="NA") { if(Enter_channel_details_now==true && marker_names_manual.length>1) //delete Hu from channel list as we are not using it for marker classification { //find index of cell_channel;; keep it as string idx_Hu=find_str_array(marker_no_manual,cell_channel); //print(idx_Hu); if(idx_Hu!="NA") //if Hu found in the channel entries, delete that corresponding channel { marker_names_manual=Array.deleteIndex(marker_names_manual, idx_Hu); marker_no_manual=Array.deleteIndex(marker_no_manual,idx_Hu); } } cell_channel=parseInt(cell_channel); if(isNaN(cell_channel)) exit("Enter channel number for cell. If leaving empty, type NA in the value"); } if(ganglia_channel!="NA") { ganglia_channel=parseInt(ganglia_channel); if(isNaN(ganglia_channel)) exit("Enter channel number for Ganglia. If leaving empty, type NA in the value"); } //Array.show(marker_names_manual); //Array.show(marker_no_manual); //exit("test"); if(sizeC>1) { if (Cell_counts_per_ganglia==true && cell_channel=="NA" && ganglia_channel=="NA") //count cells per ganglia but don't know channels for ganglia or neuron { waitForUser("Note the channels for neuron and ganglia and enter in the next box."); Dialog.create("Choose channels for "+cell_type); Dialog.addNumber("Enter "+cell_type+" channel", 3); Dialog.addNumber("Enter channel for segmenting ganglia", 2); Dialog.show(); cell_channel= Dialog.getNumber(); ganglia_channel=Dialog.getNumber(); Stack.setChannel(cell_channel); resetMinAndMax(); Stack.setChannel(ganglia_channel); resetMinAndMax(); } else if(Cell_counts_per_ganglia==true && cell_channel!="NA" && ganglia_channel=="NA") //count cells per ganglia but don't know channels for ganglia { waitForUser("Note the channels for ganglia and enter in the next box."); Dialog.create("Choose channel for ganglia"); Dialog.addNumber("Enter channel for segmenting ganglia", 2); Dialog.show(); //cell_channel= Dialog.getNumber(); ganglia_channel=Dialog.getNumber(); //Stack.setChannel(cell_channel); //resetMinAndMax(); Stack.setChannel(ganglia_channel); resetMinAndMax(); } else if(Cell_counts_per_ganglia==true && cell_channel=="NA" && ganglia_channel!="NA") //count cells per ganglia but don't know channels for neuron { waitForUser("Note the channels for "+cell_type+" and enter in the next box."); Dialog.create("Choose channel for "+cell_type); Dialog.addNumber("Enter "+cell_type+" channel", 3); Dialog.show(); cell_channel= Dialog.getNumber(); Stack.setChannel(cell_channel); resetMinAndMax(); } } //added option for extended depth of field projection for widefield images if(sizeZ>1) { print(img_name+" is a stack"); roiManager("reset"); projection_method=getBoolean("3D stack detected. Which projection method would you like?", "Maximum Intensity Projection", "Extended Depth of Field (Variance)"); if(projection_method==1) { waitForUser("Note the start and end of the stack.\nPress OK when done"); Dialog.create("Choose slices"); Dialog.addNumber("Start slice", 1); Dialog.addNumber("End slice", sizeZ); Dialog.show(); start=Dialog.getNumber(); end=Dialog.getNumber(); run("Z Project...", "start="+start+" stop="+end+" projection=[Max Intensity]"); max_projection=getTitle(); } else { max_projection=extended_depth_proj(img_name); } } else { print(img_name+" has only one slice, using as max projection"); max_projection=getTitle(); } max_save_name="MAX_"+file_name; //Segment Neurons selectWindow(max_projection); run("Select None"); run("Remove Overlay"); //if more than one channel, set on cell_channel or reference channel if(sizeC>1) { Stack.setChannel(cell_channel); } roiManager("show none"); run("Duplicate...", "title="+cell_type+"_segmentation"); seg_image=getTitle(); roiManager("reset"); n_tiles=2; new_width=round(width*scale_factor); if(new_width>1200) n_tiles=4; else if(new_width>4000) n_tiles=10; //scale image if scaling factor is not equal to 1 if(scale_factor!=1) { selectWindow(seg_image); new_width=round(width*scale_factor); new_height=round(height*scale_factor); run("Scale...", "x=- y=- width="+new_width+" height="+new_height+" interpolation=None create title=img_resize"); close(seg_image); selectWindow("img_resize"); seg_image=getTitle(); } roiManager("UseNames", "false"); selectWindow("Log"); print("*********Segmenting cells using StarDist********"); //segment neurons using StarDist model segment_cells(max_projection,seg_image,neuron_model_path,n_tiles,width,height,scale_factor,neuron_seg_lower_limit,probability,overlap); close(seg_image); //manually correct or verify if needed waitForUser("Correct "+cell_type+" ROIs if needed. You can delete or add ROIs using ROI Manager"); cell_count=roiManager("count"); rename_roi(); //rename ROIs roiManager("deselect"); selectWindow(max_projection); //uses roi to label macro code; clij is a dependency runMacro(roi_to_label); wait(5); neuron_label_image=getTitle(); //using this image to detect neuron subtypes by label overlap selectWindow(neuron_label_image); saveAs("Tiff", results_dir+"Neuron_label_"+max_save_name); rename("Neuron_label"); //saving the file will change the name, so renaming it and getting image name again neuron_label_image=getTitle(); selectWindow(neuron_label_image); run("Select None"); print("No of "+cell_type+" in "+max_projection+" : "+cell_count); roiManager("deselect"); roi_location=results_dir+cell_type+"_ROIs_"+file_name+".zip"; roiManager("save",roi_location ); selectWindow(table_name); Table.set("File name",row,file_name); Table.set("Total "+cell_type, row, cell_count); //set total count of neurons after nos analysis if nos selected Table.update; //Segment ganglia selectWindow(max_projection); run("Select None"); run("Remove Overlay"); if (Cell_counts_per_ganglia==true) { roiManager("reset"); if(Ganglia_detection=="DeepImageJ") { //ganglia_binary=ganglia_deepImageJ(max_projection,cell_channel,ganglia_channel); args=max_projection+","+cell_channel+","+ganglia_channel; //get ganglia outline runMacro(segment_ganglia,args); wait(5); ganglia_binary=getTitle(); draw_ganglia_outline(ganglia_binary,true); } else ganglia_binary=draw_ganglia_outline(ganglia_img,false); args=neuron_label_image+","+ganglia_binary; //get cell count per ganglia runMacro(ganglia_cell_count,args); //make ganglia binary image with ganglia having atleast 1 neuron selectWindow("label_overlap"); //getMinAndMax(min, max); setThreshold(1, 65535); run("Convert to Mask"); resetMinAndMax; close(ganglia_binary); selectWindow("label_overlap"); rename("ganglia_binary"); selectWindow("ganglia_binary"); ganglia_binary=getTitle(); selectWindow("cells_ganglia_count"); cell_count_per_ganglia=Table.getColumn("Cell counts"); roiManager("deselect"); ganglia_number=roiManager("count"); roi_location=results_dir+"Ganglia_ROIs_"+file_name+".zip"; roiManager("save",roi_location ); roiManager("reset"); selectWindow(table_name); Table.set("No of ganglia",0, ganglia_number); Table.setColumn("Neuron counts per ganglia", cell_count_per_ganglia); Table.update; selectWindow("cells_ganglia_count"); run("Close"); } //neuron_subtype_matrix=0; no_markers=0; //if user wants to enter markers before hand, can do that at the start //otherwise, option to enter them manually here if(marker_subtype==1) { arr=Array.getSequence(sizeC); arr=add_value_array(arr,1); if(Enter_channel_details_now==true) { channel_names=marker_names_manual;//split(marker_names_manual, ","); channel_numbers=marker_no_manual;//split(marker_no_manual, ","); channel_numbers=convert_array_int(marker_no_manual); no_markers=channel_names.length; //Array.show(channel_names); //Array.show(channel_numbers); } else { no_markers=getNumber("How many markers would you like to analyse?", 1); string=getString("Enter names of markers separated by comma (,)", "Names"); channel_names=split(string, ","); if(channel_names.length!=no_markers) exit("Channel names do not match the no of markers"); channel_numbers=newArray(sizeC); marker_label_img=newArray(sizeC); Dialog.create("Select channels for each marker"); for(i=0;i1) { channel_combinations=combinations(channel_names); //get all possible combinations and adds an underscore between name labels if multiple markers channel_combinations=sort_marker_array(channel_combinations); } else { channel_combinations=channel_names; // pass single combination //channel_combinations=channel_numbers[0]; } channel_position=newArray(); marker_label_arr=newArray(); //store names of label images generated from StarDist //count=0; //Array.show(channel_numbers); //Array.show(channel_names); //Array.show(channel_numbers); selectWindow(max_projection); Stack.setDisplayMode("color"); row=0; for(i=0;i1) marker_label_arr[i]=label_marker; marker_count=roiManager("count"); // in case any neurons added after manual verification of markers selectWindow(table_name); //Table.set("Total "+cell_type, row, cell_count); Table.set("Marker Combinations", row, channel_name); Table.set("Number of cells per marker combination", row, marker_count); //Table.set(""+cell_type, row, marker_count/cell_count); Table.update; row+=1; //selectWindow(max_projection); roiManager("deselect"); //roi_file_name= String.join(channel_arr, "_"); roi_location_marker=results_dir+channel_name+"_ROIs.zip"; roiManager("save",roi_location_marker); close(seg_marker_img); roiManager("reset"); //Array.print(marker_label_arr); if (Cell_counts_per_ganglia==true) { selectWindow(label_marker); run("Remove Overlay"); run("Select None"); args=label_marker+","+ganglia_binary; //get cell count per ganglia runMacro(ganglia_cell_count,args); close("label_overlap"); selectWindow("cells_ganglia_count"); cell_count_per_ganglia=Table.getColumn("Cell counts"); selectWindow(table_name); Table.setColumn(channel_name+" counts per ganglia", cell_count_per_ganglia); Table.update; selectWindow("cells_ganglia_count"); run("Close"); roiManager("reset"); } } //if more than one marker to analyse; if more than one marker, then it multiplies the marker labels from above to find coexpressing cells else if(channel_arr.length>=1) { for(j=0;j0) { clippedIdx = Array.trim(roiIdx,k); //roiManager("select", clippedIdx); } //else roiManager("deselect"); return k; } //add a value to every element of an array function add_value_array(arr,val) { for (i = 0; i < arr.length; i++) { temp=arr[i]+val; arr[i]=parseInt(temp); } return arr; } //convert arr to int function convert_array_int(arr) { for (i = 0; i < arr.length; i++) { arr[i]=parseInt(arr[i]); } return arr; } //function to scale imag function scale_image(img,scale_factor,name) { if(scale_factor!=1) { selectWindow(img); Stack.getDimensions(width, height, channels, slices, frames); new_width=round(width*scale_factor); new_height=round(height*scale_factor); run("Scale...", "x=- y=- width="+new_width+" height="+new_height+" interpolation=None create title="+name+"_resize"); close(img); //selectWindow(name+"_resize"); scaled_img=name+"_resize"; } else { scaled_img=img; } return scaled_img; } //generate every possible combination of markers function combinations(arr) { len=arr.length; str=""; p=Math.pow(2,len); arr_str=newArray(); for(i=0;i0) //bitwise AND comparison { //print(arr[j]); str+=arr[j]+","; } //else print("Nothing"); } if(str!="") str=substring(str,0,str.length-1); //if not empty string, remove the comma at the end arr_str[i]=str; str=""; //str+="\n"; } return arr_str; } //sort the string array based on the number of strings/markers function sort_marker_array(arr) { //print(no_combinations); rank_idx=1; rank_array=newArray(); no_markers=1; //first value is empty string, so deleting that arr=Array.deleteValue(arr,""); no_combinations=arr.length; do { for (i = 0; i < no_combinations; i++) { arr_str=split(arr[i], ","); if(arr_str.length==no_markers) { rank_array[i]=rank_idx; rank_idx+=1; } } no_markers+=1; } while (rank_idx<=no_combinations) //Array.show(arr); //Array.show(rank_array); //change order of markers based on the order specified in rank_array Array.sort(rank_array,arr); //Array.show(arr1); return arr; } //find if a string is contained within an array of strings //case insensitive function find_str_array(arr,name) { name=".*"+toLowerCase(name)+".*"; no_str=arr.length; position="NA"; for (i=0; i1) { for(ch=1;ch<=channels;ch++) { selectWindow(img); Stack.setChannel(ch); getLut(reds, greens, blues); Ext.CLIJ2_push(img); radius_x = 2.0; radius_y = 2.0; sigma = 10.0; proj_img="proj_img"+ch; Ext.CLIJ2_extendedDepthOfFocusVarianceProjection(img, proj_img, radius_x, radius_y, sigma); Ext.CLIJ2_pull(proj_img); setLut(reds, greens, blues); //Ext.CLIJ2_pull(img); concat_ch=concat_ch+"c"+ch+"="+proj_img+" "; } Ext.CLIJ2_clear(); //print(concat_ch); run("Merge Channels...", concat_ch+" create"); Stack.setDisplayMode("color"); } else { selectWindow(img); getLut(reds, greens, blues); Ext.CLIJ2_push(img); radius_x = 2.0; radius_y = 2.0; sigma = 10.0; proj_img="proj_img"; Ext.CLIJ2_extendedDepthOfFocusVarianceProjection(img, proj_img, radius_x, radius_y, sigma); Ext.CLIJ2_pull(proj_img); setLut(reds, greens, blues); } max_name="MAX_"+img; rename(max_name); close(img); return max_name; }