Demo script showing a simple method of assigning classifications for a multichannel fluorescence image

QuPath-Multiple cell classifications.groovy

/**
 * Demo script showing a simple method of assigning classifications for a multichannel fluorescence image.
 *
 * Note that at the time of writing, this requires the most recent (pre-release) QuPath version.
 * See https://petebankhead.github.io for more information.
 *
 * The channel names here refer to the sample image 'LuCa-7color_[13860,52919]_1x1component_data.tif'
 * from Perkin Elmer, after running 'Update measurement names.groovy' script.
 *
 * The original image is part of Bio-Formats demo data, under Creative Commons Attribution 4.0 International License
 * http://creativecommons.org/licenses/by/4.0/
 *
 * @author Pete Bankhead
 */

import qupath.lib.objects.PathObject
import qupath.lib.objects.classes.PathClassFactory
import qupath.lib.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 kHigh = 8
double kLow = 6
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 = [
        kThreshold('PDL1', 'Cytoplasm: PDL1 mean', kLow),
        kThreshold('CD8', 'Nucleus: CD8 mean', kHigh),
        kThreshold('FoxP3', 'Nucleus: FoxP3 mean', kHigh),
        kThreshold('CD68', 'Cytoplasm: CD68 mean', kLow),
        kThreshold('PD1', 'Cell: PD1 mean', kLow),
        kThreshold('CK', 'Cytoplasm: CK mean', kLow)
]

// 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)
        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 'None: \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)
        }

    }

}