Discussion and Script: What is a hotspot?

Hi all,
Identifying dense groups of cells, or specific microenvironment clusters is an increasingly important area of image analysis, and one that I wanted to open up a discussion on as far as automating ways of finding these types of areas.

I started with the simplest option for an earlier project, the Delaunay cluster features, which does a nice job of finding “clusters” of cells when there are only two classes of cells, or no classes. I could fairly easily find clusters of a certain size by removing all cells except for those of the class I was interested in, and then running the Delaunay analysis on each class in turn. Unfortunately, that also picked up long strings of cells, as the only requirement for a cluster was that any two cells of the cluster be within X distance of each other. In other words, it was not the measure of density I was looking for.

Aside: The reason I needed to remove every other cell class was that any two cells of a class created a line between them, a line which could not then be crossed by another class. So, if I were interested in clusters of T cells in a tumor environment, and there were more tumor cells than T-cells, the chances were I would never pick up the clusters I was interested in.
See this post about nearest neighbors for a visual demonstration of the problem: QuPath Script: Nearest Neighbors

That led to a second attempt, included here, where I built off of Pete’s script here.
https://forum.image.sc/t/find-highest-staining-region-of-a-slide/31251/7?u=mike_nelson

Script included in second post.

The new script creates a heatmap of each class, and then create isolines based on some threshold that could be turned into annotation objects.
Hopefully a clearer explanation in individual steps.

  1. Create one blank 2D array of sufficient size to represent the whole image, for each class (Negative is excluded by default, others could be as well) of interest.
  2. Cycle through all cells and add a point to the array of the appropriate class at the XY coordinates of that cell.
  3. Based on some distance, expand each of those points into a circle. If the distance is large enough, these circles will start to overlap.
  4. Create isolines at a certain value. IE, if 4 or more cells are within “distance” of a certain point, that point will be included within an annotation.
  5. Remove all annotations that contain fewer than a set number of cells. For whole tissue images, someone might be interested only in clusters of over 100 cells, while someone working with high powered confocal images might be interested in a cluster of 5-10 cells.
  6. Finish with the annotations inserted into the hierarchy so that the user can look at what other cell types are included within that annotation.

If you are only interested in one class, there is a line in the script to edit the classes that are or are not analyzed (defaults to “not Negative”). There is a smoothing value at the top of the script that can be adjusted.

Brightfield, DAB, positive cell detection.
image
image
image

Specifically looking for clusters over 30 cells, in this case. With Positive cell detection, each hotspot annotation comes with positive and negative counts, percent positive, and the area and positive cell density!
With the heatmap checkbox selected, you can also see the actual heatmap used to generate the annotations for the current run, which can help you better choose your density threshold.
image

Here is another example using the same image Pete used for his multiplex classification tutorial. With multiple classes of cells, the results become more interesting and useful.
image
image
Looking for relatively smaller clusters this time, I only found three classes with sufficient density (Teal CK, Red PDL1, Green PD1). The script currently looks for final classes, not single channel classes, so some of the areas that look like they should be within a hotspot are double or triple positive with one of the other five channels.
image

Selecting one of the red PDL1 hotspots in the top left, I can quickly see that the other main classes within the hotspot are PDL1+PD1 positive and PDL1+CD8 positive (and negative).
image
Here is the PDL1 heatmap, for reference. Each circle has a radius of 20um. That is where the “Distance between cells” measurement is used to build the map.
image
Increasing the Density to 4 removes some of the less tightly grouped hotspots and tightens the boundaries on the remaining hotspots.
image

I can use the heatmap to choose a Density value of 4 by looking at the pixel “values” under the mouse cursor on the ImageJ bar.

This script, and Pete’s that was linked earlier, are two ways of handling hotspot/cluster detection. The original script created a single circle around a hotspot, and I first expanded that script to find any locations where the density was “sufficiently high” and created circles at those points. That created a rather large number of circles, and rather than trying to fix that, I realized the circles could be built up like a topographical map (not a new idea, but new to me at the time).

From what I have seen of use cases, though, the power of being able to do this in an automated fashion relates to the meaningfulness of information that can be extracted from each hotspot. I feel that that information is more meaningful when the annotation the information is based on traces the biologically important regions.

I am interested in other ideas concerning what hotspots or clusters could mean, or other ways of finding groups of objects. This method does not capture any information about the actual structure of the microenvironment, with each class being considered separately. One possibility for capturing this type of information for a multiplex analysis might be to create isolines per marker, and use the intersection of the annotations to keep overlapping regions of certain types (T-cell -memory + Macrophage+CK) or (+memory+stroma). Another possibility with the new classification system in M9 might be to try to add sub-classes based on nearby neighbors (Tumor cells with >N T-cells within 20um) and create hotspots where the classification of the cells is based on the local environment. Combining this script with the nearest neighbors counts could open up some interesting options.

Open to ideas, hoping to get some feedback!

The script has been done for a while, but I’ve spent some time cleaning it up a bit. If anyone has ideas on how to make it run faster on large images (aside from more downsampling), I am all ears!

// Script to find high density areas of classified cells in QuPath v0.2.0-m8. Version 4.
// Expected input: Classified cells, normally within some kind of annotation that does not have the same class as the cells.
// Expected output: No change to initial cells or classifications, but the addition of classified annotations around hotspots
// Downstream: Add further measurements to annotations based on the density and percentages of classified cells?
// Script by Mike Nelson, 1/15/2020.

double smooth = 7.0

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 qupath.imagej.processing.RoiLabeling;
import ij.gui.Wand;
import java.awt.Color
import ij.measure.Calibration;
import ij.IJ
import qupath.imagej.gui.IJExtension
import qupath.lib.regions.RegionRequest
import ij.ImagePlus;
import qupath.lib.regions.ImagePlane
import qupath.lib.objects.PathObjects
import ij.process.ImageProcessor;
import qupath.lib.objects.PathCellObject
import qupath.lib.roi.ShapeSimplifier

//Default values for dialog box, or values for running the script across a project.
int pixelDensity = 3
double radiusMicrons = 20.0
int minCells = 30
boolean smoothing = true
boolean showHeatMap = false
//Collect some information from the user to use in the hotspot detection

///////////////////////////////////////////////////////////////////
def params = new ParameterList()
        
        .addIntParameter("minCells", "Minimum cell count", minCells, "cells", "Minimum number of cells in hotspot")
        //.addDoubleParameter("pixelSizeMicrons", "Pixel size", pixelSizeMicrons, GeneralTools.micrometerSymbol(), "Choose downsampling-can break script on large images if not large enough")
        .addIntParameter("density", "Density", pixelDensity, "Changes with the other variables, requires testing", "Integer values: lower pixel size requires lower density")
        .addDoubleParameter("radiusMicrons", "Distance between cells", radiusMicrons, GeneralTools.micrometerSymbol(), "Usually roughly the distance between positive cell centroids")
        .addBooleanParameter("smoothing", "Smoothing?           ", smoothing, "Do you want smoothing")
        .addBooleanParameter("heatmap", "Show Heatmap?           ", showHeatMap, "Open a new window showing the heatmap. If ImageJ is already open, you can use that to look at pixel values")
if (!Dialogs.showParameterDialog("Parameters. WARNING, can be slow on large images with many clusters", params))
    return


radiusMicrons = params.getDoubleParameterValue("radiusMicrons")
minCells = params.getIntParameterValue("minCells")
pixelDensity = params.getIntParameterValue("density")
smoothing = params.getBooleanParameterValue("smoothing")
showHeatMap = params.getBooleanParameterValue("heatmap")

///////////////////////////////////////////////////////////////////

//Comment out the entire section above and put the values you want in manually if you want to run the script "For Project"


int z = 0
int t = 0
def plane = ImagePlane.getPlane(z, t)
def imageData = getCurrentImageData()
def server = imageData.getServer()
def cells = getCellObjects()

double pixelCount = (double)server.getWidth()*(double)server.getHeight()
downsample = Math.ceil(pixelCount/(double)500000000)
pixelSizeMicrons = downsample*server.getPixelCalibration().getAveragedPixelSizeMicrons()

int w = Math.ceil(server.getWidth() / downsample)
int h = Math.ceil(server.getHeight() / downsample)
int nPixels = w * h
double radius = radiusMicrons / pixelSizeMicrons
//println("downsample " + downsample)
//println("radius " +radius)
//Unsure about this part. Maybe it shouldnt start at 0,0 but should get the upper left pixel using the imageserver?
Calibration calIJ = new Calibration();
calIJ.xOrigin = 0/downsample;
calIJ.yOrigin = 0/downsample;

//Find all classes
Set classSet = []
for (object in getCellObjects()) {
    classSet << object.getPathClass()
}
//convert classes into a list, which is ordered

/*************************************************
CLASS LIST MIGHT BE MODIFIABLE FOR MULTIPLEXING
*****************************************************/
List classList = []
classList.addAll(classSet.findAll{ 
    //If you only want one class, use it == getPathClass("MyClassHere") instead
    it != getPathClass("Negative") 
})

removeObjects(getAnnotationObjects().findAll{(classList.contains(it.getPathClass()))},true)

print("Class list: "+ classList)
println("This part may be QUITE SLOW, with no apparent sign that it is working. Please wait for the 'Done' message.")

// Create centroid map
/*****************************
Create an array of floatprocessors per class

***************************/
def fpList = []
for (aClass in classList){
    fpArray = new FloatProcessor(w,h) 
    fpList << fpArray
}
def fpNegative = new FloatProcessor(w, h)

//////////////////////// Update valid mask
//Checking for areas to ignore (outside of annotations, near borders
ByteProcessor bpValid
def annotations = getAnnotationObjects()
if (annotations) {
//making an image instead of a byteprocessor
    def imgMask = new BufferedImage(w, h, BufferedImage.TYPE_BYTE_GRAY)
        //Not sure
    def g2d = imgMask.createGraphics()
        //scale the image down by the downsample
    g2d.scale(1.0/downsample, 1.0/downsample)
        //Whole image starts black, 0s, fill annotation with white, 255?
    g2d.setColor(Color.WHITE)
    for (annotation in annotations) {
        def shape = annotation.getROI().getShape()
        g2d.fill(shape)
    }
    g2d.dispose()
        //ok, I think at this point we have a large white block defining the annotation are
    bpValid = new ByteProcessor(imgMask)
    bpValid = SimpleThresholding.thresholdAbove(new EDM().makeFloatEDM(bpValid, 0, true), radius/4 as float)
        //Ok, now we have a distance transform from the edge of the annotation object...
    //Ahah! One that is thresholded so that we don't look for hotspots near the edge of something. Not sure if I will want this behavior or not
    
}
//clear out the original annotations to make it easier to cycle through all new annotations
removeObjects(annotations,true)


///////////////////////////////////////////
//Create cell count map
for (cell in cells) {
    def roi = PathObjectTools.getROI(cell, true)
    if (roi.isEmpty())
        continue
    def pathClass = cell.getPathClass()
    
    //Ignore unclassified cells
    if (pathClass == null )
        continue
    
    int x = (roi.getCentroidX() / downsample) as int
    int y = (roi.getCentroidY() / downsample) as int
    //find the pixel of the current roi center, and validate the position against the mask that checks for being too close to the border or outside of annotations
    int check = w*y+x
    
    //This is where the fpList[] starts to get information from individual classes
    //add 1 pixel value to the fpList equivalent to the class of the cell
    //After this, each fpList object should be an array that shows COUNTS for cells within an area determined by downsampling
    
    //Make sure we are writing to the correct position in fpList
    for (i=0; i<classList.size(); i++){
        if (pathClass == classList[i]){
            if (bpValid.getf(check) != 0f){
                fpList[i].setf(x, y, fpList[i].getf(x, y) + 1f as float)
            }
        }
    }

    if (PathClassTools.isNegativeClass(pathClass) && bpValid.getf(check) != 0f)
        fpNegative.setf(x, y, fpNegative.getf(x, y) + 1f as float)

}
//At this point we have cycled through all of the cells and built N heatmaps, though they are downsampled

////////////////////////////////////////////////////////////
// In this section we create a mean filter to run across our downsampled density map, using the radius given by the user.
// This, along with the downsample, will fill in the spaces between cells
def rf = new RankFilters()

//Get an odd diameter so that there is a center
int dim = Math.ceil(radius * 2 + 5)
def fpTemp = new FloatProcessor(dim, dim)
//generate an empty square (0s) with R^2 as the center pixel value
fpTemp.setf(dim/2 as int, dim/2 as int, radius * radius as float)
//spread the radius squared across a circle using the euclidean distance, radius
rf.rank(fpTemp, radius, RankFilters.MEAN)
def pixels = fpTemp.getPixels() as float[]
//count the number of pixels within fpTemp that will actually be used by RankFilters.Mean when passing "radius"
double n = Arrays.stream(pixels).filter({f -> f > 0}).count()


// Compute sum of elements
//rankfilter is used to run a mean filter across the fpTemp area

/*######## NEED TO MAKE THIS ONLY USE INTERESTING CLASSES ###########*/
fpList.each{
    rf.rank(it, radius, RankFilters.MEAN)
    it.multiply(n)
}


//Here we take the mean-filtered density maps, apply the user's density threshold, 
for (l=0; l<fpList.size(); l++){
    //create a mask based on the user threshold
    hotspotMaskMap = SimpleThresholding.thresholdAbove(fpList[l], (float)pixelDensity)
    //not 100% sure how this line worked, but it was necessary for the getFilledPolygonROIs to function
    hotspotMaskMap.setThreshold(1, ImageProcessor.NO_THRESHOLD, ImageProcessor.NO_LUT_UPDATE)
    //use the mask to generate ROIs that surround 4 connected points (not diagonals)
    hotspotROIs = RoiLabeling.getFilledPolygonROIs(hotspotMaskMap, Wand.FOUR_CONNECTED);
    
    //print(hotspotROIs.size())

    allqupathROIs = []
    qupathROIs = []
    //convert the ImageJ ROIs to QuPath ROIs
    hotspotROIs.each{allqupathROIs << IJTools.convertToROI(it, calIJ, downsample, plane)}
    //Use the QuPath ROIs to generate annotation objects (possibly smoothed), out of the heatmap ROIs
    objects = []
    qupathROIs = allqupathROIs.findAll{it.getArea() > (radiusMicrons*radiusMicrons/(server.getPixelCalibration().getAveragedPixelSizeMicrons()*server.getPixelCalibration().getAveragedPixelSizeMicrons()))}
    smoothedROIs = []
    qupathROIs.each{smoothedROIs << ShapeSimplifier.simplifyShape(it, downsample*2)}
    
    //println("sizes "+ qupathROIs.size)
    smoothedROIs.each{objects << PathObjects.createAnnotationObject(it, classList[l]);}
 
    addObjects(objects)

}

resolveHierarchy()

//remove small hotspots
getAnnotationObjects().each{
    currentClass = it.getPathClass()
    if (classList.contains(it.getPathClass())){
        count = []
        qupath.lib.objects.PathObjectTools.getDescendantObjects(it,count, PathCellObject)
        count = count.findAll{cell -> cell.getPathClass() == currentClass}
        if (count.size < minCells){
            //print count.size
            removeObject(it,true)
        }
    }
}
Set hotSpotClassList = []
for (object in getAnnotationObjects()) {
    hotSpotClassList << object.getPathClass()
}
IJExtension.getImageJInstance()
if (showHeatMap){
    for (l=0; l<fpList.size(); l++){
        if (hotSpotClassList.contains(classList[l])){
            new ImagePlus(classList[l].toString()+" heatmap at "+ pixelSizeMicrons+ "um pixel size", fpList[l]).show()
        }
    }
}


if(smoothing){
    
    before = getAnnotationObjects()
    selectAnnotations()
    runPlugin('qupath.lib.plugins.objects.DilateAnnotationPlugin', '{"radiusMicrons": '+smooth+',  "lineCap": "Round",  "removeInterior": false,  "constrainToParent": false}');
    removeObjects(before,true)
    expanded = getAnnotationObjects()
    selectAnnotations()
    runPlugin('qupath.lib.plugins.objects.DilateAnnotationPlugin', '{"radiusMicrons": '+(-1*smooth)+',  "lineCap": "Round",  "removeInterior": false,  "constrainToParent": false}');
    removeObjects(expanded,true)
    resetSelection();
    
}
    

//return the original annotations
addObjects(annotations)
resolveHierarchy()
getAnnotationObjects().each{it.setLocked(true)}

println("Done")
1 Like

Quick update to the script posted to prevent an ROI error when downsampling on large images.
If anyone sees anything like this, please message me.