Find highest staining region of a slide

I’m trying to find the highest staining region of a WSI, but am not too sure exactly how to do it in QuPath 0.2.0. First I highlight the tissue with the wand, then run cell detection, estimate stain vector, and positive cell detection. From this point I would like to run a program that returns an area of high magnification (a rectangle with an area of 0.152 mm^2) with the highest nuclear staining. Thank you.

As a warning, you are almost certainly going to find something black. Be it a trick of the light, or a piece of dust, or some other contaminant, at some point you are going to find a very small area that is black. Though that is for brightfield, and you haven’t said what type of image you are looking at.

And just to clarify, 0.152 mm^2 is a fairly large area, about the size of a nucleus, so if you want something like that, can you not use the mean nuclear intensity/OD for your stain? Finding an exact square would probably require calculating the Square Tiles measurements for all of your cells (Analyze->Calculate features), and then sorting.
image

The average nuclear intensity is built in.

The slides I am working on look pretty clean, and I try to check for dust or other contaminants that could throw off the positive cell detection before doing anything. I am looking at brightfield H-DAB images (Ki67-stained slides).

My Compute intensity features doesn’t look the same as yours, and I can’t figure out what this feature does after running it. I used optical density in the positive cell classification also.

Making tiles that are too large misses out on some areas of the slide. Maybe it would be better to create overlapping rectangles over the larger main annotation, calculate the percentage of positive cells in those areas, and return the highest percentage rectangle.

I have the field diameter written down on my microscope, but I thought it was around 0.152 mm^2. I think I use a 22 mm eyepiece, then with a 40x objective that would be maybe even larger at 0.238 mm^2 for a single 40x field. It does look a little big when I try to make tiles that large though.

Hmm, I’m not sure I understand the end goal. Are you looking for hotspots? That was the point of the Add intensity features when you said you were detecting cells.

Min/Max aren’t terribly useful, as they only represent a single pixel. If you only want the darkest pixel on your image, then sure, you could create tiles, and add intensity features to those, though for that you just use ROI, not Square Tiles. If you want the darkest overall tile, you are going to want mean or median values.
If you want the darkest nucleus, or pixel inside of a nucleus, you are best off using the cell detection and sorting by max DAB intensity. If you want to get a bit more exact and stain free, go with max OD sum from Add Intensity Features (if you haven’t calibrated your DAB in the Preprocessing->Estimate stain vectors). You can sort in the Measure->Show detection measurements. In that case you should also go with ROI.

In this image, the selected cell has the darkest pixel.


That cell also has XY coordinates you could use to do other things with, like Objects->Annotation->Specify annotation
The measurements will depend on your pixel size.

If you use M5, you should have an additional option to include only the Nucleus as your region of interest in the same dropdown with Square Tiles, etc. when using Add Intensity Features.

Yes, I am looking for hotspots. I want to find an area of a 40x magnification field that has the highest staining percentage within a larger annotation that I made with the wand.
This is what I’m starting with:

And I would like to find the hotspot with the highest percentage of staining and have that be an annotation, something like this:

Eyeballing it and using the specify annotation tool to make a rectangle may work for now, but I was just wondering if there was an easy way to compute this.

I’m still not sure how the Compute intensity feature works, doesn’t seem to do anything after I run it on any of the settings.

It depends on what objects you are running the intensity features on. It adds those features to the objects, and is a general tool for classification. More information in the guide here.

There is already a script that does pretty much exactly what you want, though I’ll need to hunt it down.

Here’s a (slightly-involved) script that I once wrote, adapted to run under QuPath v0.2.0-m5.

// Script to find highest density of positive cells in QuPath v0.2.0-m5.
// Assumes cell detection has already been calculated.
// Algorithm runs at a specified resolution (here 20 um/px).

import ij.plugin.filter.EDM
import ij.plugin.filter.RankFilters
import ij.process.Blitter
import ij.process.ByteProcessor
import ij.process.FloatProcessor
import qupath.imagej.processing.SimpleThresholding
import qupath.lib.objects.classes.PathClass
import qupath.lib.objects.classes.PathClassFactory
import qupath.lib.plugins.parameters.ParameterList

import java.awt.Color

double pixelSizeMicrons = 20.0
double radiusMicrons = 1000.0
int minCells = 500
boolean tumorOnly = false


def params = new ParameterList()
        .addDoubleParameter("radiusMicrons", "Radius", radiusMicrons, GeneralTools.micrometerSymbol(), "Choose radius of hotspot")
        .addIntParameter("minCells", "Minimum cell count", minCells, "", "Minimum number of cells in hotspot")
        .addBooleanParameter("limitTumor", "Limit to tumor cells", tumorOnly, "Only consider tumor cells in density calculations")

if (!Dialogs.showParameterDialog("Parameters", params))
    return

radiusMicrons = params.getDoubleParameterValue("radiusMicrons")
minCells = params.getIntParameterValue("minCells")
tumorOnly = params.getBooleanParameterValue("limitTumor")


def filterClass = tumorOnly ? getPathClass('Tumor') : null

def imageData = getCurrentImageData()
def server = imageData.getServer()
def cells = getCellObjects()

double downsample = pixelSizeMicrons / server.getPixelCalibration().getAveragedPixelSizeMicrons()
int w = Math.ceil(server.getWidth() / downsample)
int h = Math.ceil(server.getHeight() / downsample)


// Create centroid map
def fpNegative = new FloatProcessor(w, h)
def fpPositive = new FloatProcessor(w, h)
def unknownClasses = new HashSet<PathClass>()
for (cell in cells) {
    def roi = PathObjectTools.getROI(cell, true)
    if (roi.isEmpty())
        continue
    def pathClass = cell.getPathClass()
    if (pathClass == null || (filterClass != null && !pathClass.isDerivedFrom(filterClass)))
        continue
    int x = (roi.getCentroidX() / downsample) as int
    int y = (roi.getCentroidY() / downsample) as int
    if (PathClassTools.isPositiveClass(pathClass))
        fpPositive.setf(x, y, fpPositive.getf(x, y) + 1f as float)
    else if (PathClassTools.isNegativeClass(pathClass))
        fpNegative.setf(x, y, fpNegative.getf(x, y) + 1f as float)
    else
        unknownClasses.add(pathClass)
}
if (!unknownClasses.isEmpty())
    print('Unknown classes: ' + unknownClasses)

// Compute sum of filter elements
def rf = new RankFilters()
double radius = radiusMicrons / pixelSizeMicrons
int dim = Math.ceil(radius * 2 + 5)
def fpTemp = new FloatProcessor(dim, dim)
fpTemp.setf(dim/2 as int, dim/2 as int, radius * radius as float)
rf.rank(fpTemp, radius, RankFilters.MEAN)
def pixels = fpTemp.getPixels() as float[]
double n = Arrays.stream(pixels).filter({f -> f > 0}).count()

// Compute sums
rf.rank(fpPositive, radius, RankFilters.MEAN)
fpPositive.multiply(n)
rf.rank(fpNegative, radius, RankFilters.MEAN)
fpNegative.multiply(n)

// Update valid mask
ByteProcessor bpValid
def annotations = getAnnotationObjects()
if (annotations) {
    def imgMask = new BufferedImage(w, h, BufferedImage.TYPE_BYTE_GRAY)
    def g2d = imgMask.createGraphics()
    g2d.scale(1.0/downsample, 1.0/downsample)
    g2d.setColor(Color.WHITE)
    for (annotation in annotations) {
        def shape = annotation.getROI().getShape()
        g2d.fill(shape)
    }
    g2d.dispose()
    bpValid = new ByteProcessor(imgMask)
    bpValid = SimpleThresholding.thresholdAbove(new EDM().makeFloatEDM(bpValid, 0, true), radius as float)
}

// Compute local densities
def fpDensity = new FloatProcessor(w, h)
def fpDensityAll = new FloatProcessor(w, h)
def fpSum = fpNegative.duplicate()
fpSum.copyBits(fpPositive,0, 0, Blitter.ADD)
int nPixels = w * h
float maxDensity = Float.NEGATIVE_INFINITY
int maxInd = -1
for (int i = 0; i < nPixels; i++) {
    float total = fpSum.getf(i)
    if (total == 0f)
        continue
    if (bpValid.getf(i) == 0f)
        continue
    float density = fpPositive.getf(i) / total
    fpDensityAll.setf(i, density)
    if (total >= minCells) {
        fpDensity.setf(i, density)
        if (density > maxDensity) {
            maxDensity = density
            maxInd = i
        }
    }
}
if (maxInd < 0) {
    Dialogs.showWarningNotification("Hotspot detection", "No hotspot found!")
    println 'No region found!'
    return
}
double x = downsample * (maxInd % w)
double y = downsample * (maxInd / w)

double fullRadius = radiusMicrons / server.getPixelCalibration().getAveragedPixelSizeMicrons()
def roi = ROIs.createEllipseROI(x - fullRadius, y - fullRadius, fullRadius*2, fullRadius*2, ImagePlane.getDefaultPlane())
def hotspot = PathObjects.createAnnotationObject(roi)

hotspot.setName("Hotspot")
addObject(hotspot)

//new ImagePlus("Density", fpDensity).show()

Basically it attempts to find the ‘hottest hotspot’ defined as the circle of a specified size with the highest density of positive cells (optionally restricted to tumor cells only).

Note that the result is not exact, because it involves creating a binned version of the image to determine the local counts efficiently.

3 Likes

That is exactly what I was looking for! It didn’t work on M4 but did on M5. The window that lets you choose the area and limit it to tumor cells is very helpful.

I am not exactly sure what you mean by the binning limitation. In which situations could this cause erroneous results?

1 Like

The default settings are reducing the pixel size (resolution) to 20um per pixel in the image that is analyzed. It shouldn’t miss any cells, but might not cluster them in a way you would expect depending on where the border of the 20um “pixel” happens to fall.

2 Likes

I mean that it does the calculations by effectively creating a new, lower-resolution image with pixels that are a certain (very large) size, e.g. 20 um, and then counting up the cells that fall within each of these pixels for a predefined area.

Whenever the result is mapped back to the full image, the actual circle is drawn at the full resolution, without this binning involved. This means that some cells might have contributed to the density calculation but fall outside the actual circle that is drawn (because they are just over the boundary… but close enough to be in the same bin as cells inside the boundary).

You can see this sometimes if you set your minimum cell count for the density calculation to be (for example) 1000, but the actual number of cells QuPath counts inside the circle is a bit less.

You can reduce this ‘error’ by adjusting the double pixelSizeMicrons = 20.0 line in the script, but choosing smaller values will mean the script takes longer to run… and probably still finds the same region (perhaps shifted by a few pixels).

I have no reason to think the results will be erroneous, only that they won’t be exact; there might be a ‘hotter hotspot’ somewhere, but it is likely to be very close to the one that was identified… and the script is still likely to be rather a lot more reliable than estimating it by eye.

(However if you do find any errors, please let me know!)

1 Like