Multichannel fluorescence detection improvement

Hi all. I am currently using QuPath to determine cell phenotypes in a 4 channel fluorescent image. I am using the DAPI channel to detect nuclei and then using this script to perform classifications (this is the same script Pete posted on his blog in this post, adapted for my phenotypes and with some color classifications added at the end).

/**
 * Script showing a simple method of assigning multiple classifications for a multichannel fluorescence image.
 *
 * Written on QuPath (0.2.0-m3)
 *
 * Note: this script should be run after an initial cell detection has been run
 *
 * @author Pete Bankhead
 */
 
import qupath.lib.objects.PathObject
import qupath.lib.objects.classes.PathClassFactory
import qupath.lib.gui.scripting.QPEx


// Check if the version is ok
print 'Warning! This script requires the latest QuPath updates (currently v0.1.3 beta)'
def version = qupath.lib.gui.QuPathGUI.getInstance().versionString
if (version == null) {
    print 'Could not find the version string, so will give it the benefit of the doubt...'
} else if (version <= '0.1.2') {
    print 'Sorry, this QuPath version is not supported!'
    return
}
print 'here'

// Create some simple classifiers based on a single measurement value,
// using either absolute thresholds or one computed using the median absolute deviation
double k_1 = 800
double k_2 = 800
double k_3 = 800
def kThreshold = {def a, def b, def c -> return MeasurementClassifier.createFromK(a, b, c)}
def absThreshold = {def a, def b, def c -> return MeasurementClassifier.createFromAbsoluteThreshold(a, b, c)}
def classifiers = [
        absThreshold('CCR6+', 'Cell: Channel 1 mean', k_1),
        absThreshold('CD3+', 'Cell: Channel 3 mean', k_3),
        absThreshold('ENV+', 'Cell: Channel 2 mean', k_2),
]

// Lets set the default color to be the same as the default corresponding channel, if we can
try {
    def server = QPEx.getCurrentImageData().getServer()
    def previousColors = []
    for (int c = 0; c < server.nChannels(); c++) {
        def channelName = server.getChannelName(c)23
        def color = server.getDefaultChannelColor(c)
        // Don't reuse colors!
        if (color in previousColors)
            continue
        for (classifier in classifiers) {
            if (channelName.contains(classifier.classificationName)) {
                classifier.defaultColor = color
                previousColors.add(color)
                break
            }
        }
    }
} catch (Exception e) {
    println 'Unable to set colors from channels: ' + e.getLocalizedMessage()
}

// Apply classifications
def cells = QPEx.getCellObjects()
cells.each {it.setPathClass(null)}
for (classifier in classifiers) {
    println classifier.classifyObjects(cells)
}
println 'Uninfected Cell: \t' + cells.count {it.getPathClass() == null}
QPEx.fireHierarchyUpdate()


/**
 * Helper class to calculate & apply thresholds, resulting in object classifications being set.
 */
class MeasurementClassifier {

    String classificationName
    String measurementName
    double threshold = Double.NaN
    double k
    Integer defaultColor

    /**
     * Create a classifier using a calculated threshold value applied to a single measurement.
     * The actual threshold is derived from the measurement of a collection of objects
     * as median + k x sigma, where sigma is a standard deviation estimate derived from the median
     * absolute deviation.
     *
     * @param classificationName
     * @param measurementName
     * @param threshold
     * @return
     */
    static MeasurementClassifier createFromK(String classificationName, String measurementName, double k) {
        def mc = new MeasurementClassifier()
        mc.classificationName = classificationName
        mc.measurementName = measurementName
        mc.k = k
        return mc
    }

    /**
     * Create a classifier using a specific absolute threshold value applied to a single measurement.
     *
     * @param classificationName
     * @param measurementName
     * @param threshold
     * @return
     */
    static MeasurementClassifier createFromAbsoluteThreshold(String classificationName, String measurementName, double threshold) {
        def mc = new MeasurementClassifier()
        mc.classificationName = classificationName
        mc.measurementName = measurementName
        mc.threshold = threshold
        return mc
    }

    /**
     * Calculate threshold for measurementName as median + k x sigma,
     * where sigma is a standard deviation estimate
     * derived from the median absolute deviation.
     *
     * @param pathObjects
     * @return
     */
    double calculateThresholdFromK(Collection<PathObject> pathObjects) {
        // Create an array of all non-NaN values
        def allMeasurements = pathObjects.stream()
                .mapToDouble({p -> p.getMeasurementList().getMeasurementValue(measurementName)})
                .filter({d -> !Double.isNaN(d)})
                .toArray()
        double median = getMedian(allMeasurements)

        // Subtract median & get absolute value
        def absMedianSubtracted = Arrays.stream(allMeasurements).map({d -> Math.abs(d - median)}).toArray()
        Arrays.sort(absMedianSubtracted)

        // Compute median absolute deviation & convert to standard deviation approximation
        double medianAbsoluteDeviation = getMedian(absMedianSubtracted)
        double sigma = medianAbsoluteDeviation / 0.6745

        // Return threshold
        return median + k * sigma
    }

    /**
     * Get median value from array (this will sort the array!)
     */
    double getMedian(double[] vals) {
        if (vals.length == 0)
            return Double.NaN
        Arrays.sort(vals)
        if (vals.length % 2 == 1)
            return vals[(int)(vals.length / 2)]
        else
            return (vals[(int)(vals.length / 2)-1] + vals[(int)(vals.length / 2)]) / 2.0
    }

    /**
     * Classify objects using the threshold defined here, calculated if necessary.
     *
     * @param pathObjects
     * @return
     */
    ClassificationResults classifyObjects(Collection<PathObject> pathObjects) {
        double myThreshold = Double.isNaN(threshold) ? calculateThresholdFromK(pathObjects) : threshold
        def positive = pathObjects.findAll {it.getMeasurementList().getMeasurementValue(measurementName) > myThreshold}
        positive.each {
            def currentClass = it.getPathClass()
            def pathClass
            // Create a classifier - only assign the color if this is a single classification
            if (currentClass == null) {
                pathClass = PathClassFactory.getPathClass(classificationName, defaultColor)
                if (defaultColor != null)
                    pathClass.setColor(defaultColor)
            }
            else
                pathClass = PathClassFactory.getDerivedPathClass(currentClass, classificationName, null)
            it.setPathClass(pathClass)
        }
        return new ClassificationResults(actualThreshold: myThreshold, nPositive: positive.size())
    }


    class ClassificationResults {

        double actualThreshold
        int nPositive

        String toString() {
            String.format('%s: \t %d \t (%s > %.3f)', classificationName, nPositive, measurementName, actualThreshold)
        }

    }

}

// Access phenotype sub-classifications
ccr6 = getPathClass('CCR6+')
env = getPathClass('ENV+')
cd3 = getPathClass('CD3+')


// Set the color, using a packed RGB value
color = getColorRGB(0, 0, 255)
ccr6.setColor(color)

// Update the GUI
fireHierarchyUpdate()

// Set the color, using a packed RGB value
color = getColorRGB(255, 0, 0)
env.setColor(color)

// Update the GUI
fireHierarchyUpdate()

// Set the color, using a packed RGB value
color = getColorRGB(0, 255, 0)
cd3.setColor(color)

// Update the GUI
fireHierarchyUpdate()

The threshold function in this script was good for initial analysis, but ultimately I want to have a more robust way to detect markers. Using guidelines set out in a paper that performed similar image analysis using different software, I would like cells to be considered positive for specific markers only if all of the following conditions are met;

  1. The mean cell intensity is greater than the upper quartile of mean cell intensity values in a specific channel (will act as a threshold).
  2. The standard deviation cell intensity if greater than the average standard deviation cell intensity in a specific channel (will solve some autofluorescence / background issues).
  3. (Optional) The mean cell intensity is lower than a specific value (will also help to solve background issues).

I have been trying to script this myself, but unfortunately I do not know groovy and am having some difficulty finding documentation on methods / functions online (along these lines, any suggestions on how to learn groovy / finding documentation would also be appreciated).

Any help on these matters would be greatly appreciated, thanks!

I don’t have time to dig into this quite yet, but aside from definitely needing some scripting, there are some warning flags, so I hope you have a good grasp on them!

  1. Imagine what this would do for a negative control slide, if it is picking this value from a single slide’s data. You would need some very specific reasoning and sample types to use this sort of moving target threshold, as all negative samples would have a lower threshold and thus be less negative.
  2. Do you want the standard deviation of the whole image here, or are you sure you only want rare events? This sounds like it would eliminate much of the positive signal in a highly variable sample.
  3. I would generally consider this to be the default threshold, with the others cherry picking values from it.

Anyway, you may have already considered all of these things, and I don’t know what your experiment is about, but do keep them in mind. Getting an accurate threshold on a moving target is… difficult when the distribution changes as well.

Thanks for the quick response!

So far, I have been exploring QuPath as a possible alternative to the current ways our lab performs cell counting (I’ve only been using the software for about a week or two), so a lot of this is just initial / test analysis to see if we can improve our detection / replicate data we have already produced via more arduous methods. In response to your points:

  1. The upper quartile intensity values and standard deviation intensity values within detected cells for each image will be scaled by a constant (which will be determined based on visualize inspection). I should have mentioned this previously. According to the paper which I am basing this analysis on, scaling these parameters by the same constant value for each image worked as well as scaling them by a unique one for each image.
  2. Standard deviation of intensity within detected cells.

I appreciate your concerns raised, and it is key to consider these issues moving forward. Your last sentence is of course the main rub; I’m not sure if this way of analysis will be the best way to go, but I thought it would be a little more dynamic / robust than cherry picking intensity values from each image to set as thresholds.

Yep, sorry I wasn’t providing anything more useful. As long as you can describe what you want in math terms (which it seems you have done), it should be able to be coded in the section where the thresholds are defined. You would essentially want a new function after:

def kThreshold = {def a, def b, def c -> return MeasurementClassifier.createFromK(a, b, c)}
def absThreshold = {def a, def b, def c -> return MeasurementClassifier.createFromAbsoluteThreshold(a, b, c)}

Something like def myThreshold, and then a function in the same place as double calculateThresholdFromK about halfway through the script. *Edit actually it looks like you would want a createFromMyThreshold function first (around 101 in Pete’s script), and then a calculate function. So two additional functions for your classifier.

You can see both mean and standard deviation calculated at various points throughout the script. Upper quartile might be trickier, as I don’t know the function to do that off the top of my head. Maybe @petebankhead will have something quick, but you may have to calculate that on your own based on sorting then getting the size of the array times 0.75-end.

I prefer the staining to be as consistent as possible, as using machine learning or calculating thresholds requires far better sample preparation than manual determinations, especially without a consistent control. Patient samples can be tricky though, and I have seen groups laboriously determine manual thresholds individually for many patient samples in order to achieve as high an accuracy as possible (hopefully blinded, and based on “that looks like CD8 staining” rather than "this patient should have high CD8 counts). In the end, accuracy is the goal.

Thank you for the advice on staining. With regards to the code, I understand the premise, but am confused about defining (in Pete’s original script) kHigh and kLow, as they seem necessary to run the function. However, I am using multiple parameters to classify as infected cells, not just one threshold value.

Maybe I am interpreting this wrong (as I stated in my initial post, I am not that well versed in groovy).

There are various ways of doing this, but I was recommending adding to the script rather than editing the kNN classifier. Essentially, add a third classifier that does whatever you want it to along with the kNN and simple fixed threshold classifier. The steps to do that based on the layout of this script involve a function that returns whether it is a particular class, and another that determines which of those classes an object is positive for. Most of your code/algorithm would go into replacing that second function, which would not be a simple single step.