IJ skeleton/foci analysis

soapes.ijm

// search for TODO sections to see parameters that still need tuning

//This macro takes a 3-channel 16 bit immunofluorescance image as input.
//  channel 1 (red): SC
//  channel 2 (green): foci
// channel 3 (blue): centromeres

//The primary goal of this macro is to identify rois covering main cell features,
//measure the total distance of skeletonized SC, number of foci on each SC and distance from foci to the centromere on the same skeleton
// Requirements:
//		Analyze>Set Measurements... includes Centroid measurement
//	 	Enable the "Biomedgroup" update site (see http://wiki.imagej.net/How_to_follow_a_3rd_party_update_site)
//Output is printed to log and results window.

// *** Start of macro ***
//centromere area ~30

// split images
run("Set Scale...", "distance=0 known=0 global");//remove previous scales
T = getTitle;
selectWindow(T);
run("Stack to Images");

// cache image names
selectImage(1);
redImage = getTitle();
selectImage(2);
greenImage = getTitle();
selectImage(3);
blueImage = getTitle();

imageCalculator("and", greenImage, redImage); //limit the foci channel to those overlapping SC
selectWindow(blueImage);

// process centromeres (blue channel)
// TODO may want to adjust threshold/dilation. The size of the centromers affects
// the accuracy in assigning them to skeletons
// minimum threshold is good to restrict area to JUST known centromeres
setAutoThreshold("Minimum dark");
run("Convert to Mask");
run("Despeckle");
// Minimum threshold is fairly aggressive so dilate once to increase centromere area.
// This will help match SCs with centromeres successfully.
run("Dilate");
run("Analyze Particles...", "  show=Outlines display exclude summarize add");
centromereX = newArray();
centromereY = newArray();
for (i=0;i<roiManager("count");i++){
	x = getResult("X", i);
	centromereX = Array.concat(centromereX, x);
	y = getResult("Y", i);
	centromereY = Array.concat(centromereX, y);
	centAi = getResult("Area", i);
	print("the area of centromere "+i+" is "+centAi);
}

setBatchMode("hide"); // hide the UI for this computation to avoid unnecessary overhead
centromereCount = roiManager("count");//centromere indices
print("number of centromeres: "+centromereCount);
for(i=0;i<roiManager("count");i++){
	roiManager("select",i);
	cIndex = i+1;
	roiManager("Rename", "centromere "+ cIndex);//renaming to keep track of rois
	setResult("Centromere index", i, cIndex); // link to centromere roi index
}

//process SC (red channel)
selectWindow(redImage);
setAutoThreshold("Default dark");
run("Convert to Mask");
run("Dilate");
run("Despeckle");
run("Skeletonize");
run("Invert LUT");

// Requires the Ridge Detection plugin from the Biomedgroup update site.
run("Ridge Detection", "line_width=2 high_contrast=255 low_contrast=240 estimate_width extend_line show_junction_points show_ids displayresults add_to_manager method_for_overlap_resolution=SLOPE sigma=1.2 lower_threshold=16.83 upper_threshold=40");
// NOTE: if you find ridge detection is not detecting your lines nicely, you can comment out the above line, and uncomment the following "waitForRser" line. This will
// allow you to properly segment your lines.
//waitForUser("Wait for Ridge Detection...", "Please run the Ridge Detection plugin.\nTune the parameters as needed - preview mode is recommended.\nWhen complete, click \"OK\" in this dialog to continue.");

setBatchMode("hide"); // hide the UI for this computation to avoid unnecessary overhead

// Clean up from Ridge Detection
// Remove junction points
scCount = 0;
for(i=centromereCount;i<roiManager("count");i++){
	roiManager("select",i);
	if (startsWith(Roi.getName(), "JP-")) {
		roiManager("delete");
		i--;
	}
	else {
		scCount++;
		roiManager("Rename", "SC "+ scCount);//renaming to keep track of rois
		roiManager("measure");
		setResult("SC index", i, scCount); // link to SC roi index
	}
}
print("number of SC skeletons: " + scCount);
setBatchMode("exit and display");

// determine which centromers are contained in SCs
for(i=0;i<centromereCount;i++) {
	found = false;

	//NB: I tried iterating over the roi bounds once to find 255 pixels and using the position of that pixel to find containing skeletons... with poor results!
	//     may be worth further investigation?
	// find the skeleton containing this centromere
	for (j=centromereCount;j<centromereCount+scCount && !found;j++) {
		roiManager("select", j); // select next SC

		pair = getResult("Paired Centromere", j);

		// skip previously paired SCs
		if (isNaN(pair) || pair==0) {
			// get coordinates of this skeleton
			Roi.getCoordinates(xPoints, yPoints);

			// select the current centromere
			roiManager("select", i);

			// Limit our search area on the skeleton.
			// We won't go more than "maxDist" away from either end point.
			Roi.getBounds(cX, cY, cW, cH);
			//TODO decide what this value should be..
			centromereTolerance = 1.7; // started with 1.7 centromere lengths
			maxDist = ((cW + cH) / 2) * centromereTolerance;
			index = -1;

			// search from index = 0;
			for (k=0; k<xPoints.length && k<maxDist && !found; k++) {
				x = xPoints[k];
				y = yPoints[k];
				if (Roi.contains(x, y)) {
					found = true;
					index = k;
				}
			}

			// search from index=xPoints.length-1
			for (k=xPoints.length-1; k>0 && k>=xPoints.length - maxDist && !found; k--) {
				x = xPoints[k];
				y = yPoints[k];
				if (Roi.contains(x, y)) {
					found = true;
					index = k;
				}
			}

			if (found) {
				setResult("Paired Centromere", j, i+1); // link SC to centromere index
				setResult("Index in skeleton", i, index); // record the index on the skeleton that the centromere was encountered.
				print("Pairing centromere " + (i+1) + " with SC " + (j+1-centromereCount));
			}
		}
	}
}

setBatchMode("exit and display");

//process foci (green channel)
selectWindow(greenImage);
//TODO play with threshold values (Image>Adjust>Threshold...) and algorithm
//TODO play with analyze particles min size
run("Auto Threshold", "method=Intermodes white");
run("Convert to Mask");
// Foci are larger than the background particles so we want to exclude small objects
run("Analyze Particles...", "size=6-Infinity display exclude summarize add");

print("Found " + (roiManager("count") - scCount - centromereCount) + " foci");

setBatchMode("hide"); // hide the UI for this computation to avoid unnecessary overhead

// label foci
fociCount = 0;
for(i=centromereCount + scCount;i<roiManager("count");i++) {
	fociCount++;
	roiManager("select", i);
	roiManager("Rename", "foci "+ fociCount);//renaming to keep track of rois
	setResult("Foci index",i, fociCount);
}

// For each skeleton, find the foci contained within
for (j=centromereCount;j<centromereCount+scCount;j++) {
	roiManager("select", j); // select next SC
	print("Looking for foci in skeleton " + (j+1-centromereCount) + " of " + scCount + "...");

	// get points for this SC
	Roi.getCoordinates(xPoints, yPoints);

	// For each point in this skeleton, search the foci for a match.
	for (k=0; k<xPoints.length; k++) {
		x = xPoints[k];
		y = yPoints[k];

		for(i=centromereCount + scCount;i<roiManager("count");i++) {
			result = getResult("Paired Centromere", i);
			// only check foci not already matched
			if (isNaN(result) || result ==0) {
				// make the current foci active
				roiManager("select", i);
				if (Roi.contains(x, y)) {
					setResult("Paired Skeleton", i, j+1-centromereCount); // record the skeleton containing this foci

					//add 1 to the spot count of the SC
					count = getResult("Foci count", j);
					if (isNaN(count)) count = 0;
					count++;
					setResult("Foci count", j, count);

					//measure euclidian distance to the centromere of that SC
					centIndex = getResult("Paired Centromere", j);
					if (centIndex > 0) {
						setResult("Paired Centromere", i, centIndex);
						setResult("Index in skeleton", i, k);
						centIndex--; // adjust for centromere/row index offset
						centPos = getResult("Index in skeleton", centIndex);

						// Set the skeletal length, which is simply the difference in
						// indices within the skeleton between centromere and foci.
						if (!isNaN(centPos)) {
							length = abs(centPos - k);
							setResult("Foci distance", i, length);
						}
					}
					print("Found foci " + (i+1) + " in SC " + (j+1-centromereCount));
					wait(50); // wait briefly when we find a match to ensure the UI has caught up.
					              // not 100% sure this is needed just had some cases of foci getting multiple matches.
				}
			}
		}
	}
}
setBatchMode("exit and display");

print("*** Centromere analysis complete");