IF Colocalization Calculations for QuPath

Hi all, another script. This one unfortunately does not have a GUI at the moment, but I wanted to get the option for colocalization out there as soon as I had it working. If you find that it is not working for you, please do let me know! It is certainly a bit rough around the edges, and doesn’t check for things like attempting to run the script on brightfield images. It is also fairly slow as it only runs in one thread. The code can be found here, and is intended for 0.2.0m2 (EDIT: Test version for 0.1.2 and 0.1.3 here. Thanks Pete for all of the help!).
How To:
First, have an image with some detections. Currently the script supports general detection objects (all), cells, subcellular objects, tiles (including SLICs), nuclei and cytoplasms. Once you have some objects, go into the script and edit several values so that the script will do what you want it to. You should see something like this:

//Channel numbers are based on cell measurement/channel order in Brightness/contrast menu, starting with 1

int FIRST_CHANNEL = 2

int SECOND_CHANNEL = 4

//CHOOSE ONE: "cell", "nucleus", "cytoplasm", "tile", "detection", "subcell"

//"detection" should be the equivalent of everything

String objectType = "cell"

//These should be figured out for a given sample to eliminate background signal

//Pixels below this value will not be considered for a given channel.

//Used for Manders coefficients only.

ch1Background = 4000

ch2Background = 1500

Pick the two channels you want to perform the colocalization analysis on; these numbers should be the same as the “Channel #” in the measurements list, or can be determined by the order in the Brightness/Display menu.

Next, select the detection type you want to analyze. “detection” would be the most general, and I think will hit everything that isn’t an annotation in the image. Tiles should get both tile objects and any kind of SLIC or other superpixel. Cells, nuclei, and cytoplasm should be fairly self explanatory, though your results might be odd if you are using one of the non-standard cell detections. Check the object types in the measurement list or ask here if you are getting odd results.

To get meaningful Manders coefficients you will need background values for each channel. The simplest way to think about these are that any pixel in a given channel below that value will be ignored when checking for overlap. I highly recommend reading up on colocalization before attempting to interpret the meaning of your results. If anyone has recommendations for short, clear, and to the point references to add here, please let me know and I will include links.
ImageJ Coloc

An example, with pictures!
image
Starting here, I have cells with subcellular detections created from channel 2. I was curious how channels 5 and 6 would colocalize, both within the channel 2 areas, and in general within the cells. All three channels are for cytoplasmic proteins.


The Pearson’s coefficient shows that the 5 and 6 markers are negatively correlated within the channel 2 areas, while the Manders cofficients show that most of the channel 6 stain shows up within the channel 5 stain… within the channel 2 subcellular regions. But is the fact that I looked in the subcellular regions important?

Looking at the cytoplasmic correlation coefficients, we can see that in general, the 5/6 markers are not strongly correlated, though the areas that were “unstained” (background only) do show high correlation. That is something to watch out for! The Manders coefficients for the cytoplasm also show that the channel 6 positive areas are subsets of the channel 5 positive areas, regardless of whether we look in the channel 2 areas alone. It would have been incorrect to draw any strong conclusion that this relationship was unique to the channel 2 subcells! It also could be that my threshold was a bit low on the channel 6 Manders cutoff. I would always recommend testing out different values, or better yet, having a negative control to set your thresholds!

Feedback/criticism/suggestions for improvements encouraged!

Update for 0.2.0m7 here.
https://gist.github.com/Svidro/68dd668af64ad91b2f76022015dd8a45#file-colocalization-of-channels-per-detection-0-2-0m7-groovy

I should note I have been a bit busy and have not had time to fully test/vet this. Please make some noise if there are any problems.

1 Like

I’ve tested this in the m9 milestone and it’s been working well for me so far. I was interested in applying the Costes method here however as a way to automatically determine the threshold. I’ve included the code below and it should be inserted right before the Manders calculations. It is also important to note that Costes method may not always be a reliable method to determine thresholds as outlined in the reference provided by Research Associate (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3074624/).

Be advised that I have not tested this thoroughly yet and am not an expert on colocalization at all. I’ve written the following based on how I interpreted their description of the algorithm and then implemented it into the code above. I may have misunderstood some parts so please correct me if any of this looks wrong!

//Costes Method
    
//Calculate the slope and intercept for the line of regression
//where slope = r(Sy/Sx) and y-intercept = Ymean - m*Xmean we get the following after simplification
    
double slope = pearsonTop / botCh1.sum()
double intercept = ch2Mean - slope*ch1Mean
    
//Incriment for lowering threshold
double incriment = 0.10

//Determine the highest pixel value for the first channel to set as the first threshold    
ch1Max = ch1Pixels[0]
for(int i = 0; i < ch1Pixels.size(); i++){
    if(ch1Pixels[i] > ch1Max) ch1Max = ch1Pixels[i]
}
    
double threshold1 = ch1Max + incriment
double threshold2

//Not sure if it is neccesary to check for negative thresholds but I will anyways
while(pearson > 0 && threshold1 > 0){
    
    //Lower thrshold
    threshold1 -= incriment
    threshold2 = slope*threshold1 + intercept
    
    //Redo all Pearson's calculations here but only using pixels below both thresholds
        
    ch1 = []
    ch2 = []
    for (i=0; i<size; i++){
        if(bpSLICs[i] && ch1Pixels[i] < threshold1 && ch2Pixels[i] < threshold2){
            ch1<<ch1Pixels[i]
            ch2<<ch2Pixels[i]
        }
    }
        
    //get the new number of pixels to be analyzed
    size = ch1.size()
    if(size==0) break
    //Calculating the mean for Pearson's
    ch1Mean = ch1.sum()/size
    ch2Mean = ch2.sum()/size
        
    //Create the sum for the top half of the pearson's correlation coefficient
    top = []
    for (i=0; i<size;i++){top << (ch1[i]-ch1Mean)*(ch2[i]-ch2Mean)}
    pearsonTop = top.sum()
    
    //Sums for the two bottom parts
    botCh1 = []
    for (i=0; i<size;i++){botCh1<< (ch1[i]-ch1Mean)*(ch1[i]-ch1Mean)}
    rootCh1 = Math.sqrt(botCh1.sum())
    
    botCh2 = []
    for (i=0; i<size;i++){botCh2 << (ch2[i]-ch2Mean)*(ch2[i]-ch2Mean)}
    rootCh2 = Math.sqrt(botCh2.sum())
        
    pearsonBot = rootCh2*rootCh1
        
    pearson = pearsonTop/pearsonBot
        
}
    
ch1Background = threshold1
ch2Background = threshold2
1 Like

I’m also not sure if the background intensity (or threshold) should be determined for each detection or set the same for every detection? In the script above, the threshold is set individually for each detection based on their individual PCC measurements so different detections will have diffferent thresholds which seems off to me. If anyone has any suggestions as to which method for determining background would be best, please let me know as I am quite new to all this and most literature I have found does not thoroughly explain how background intensity was determined, despite MCC being extremely sensitive to background intensity settings.

Oh, right, it’s generally used to determine the overall image threshold. Which is part of why I’m usually not a fan, but you should run it using every pixel in the image. The problem is that might require an array that is too large for Java to handle, which is why QuPath is so handy vs ImageJ in general (in terms of being able to open such images). If your images are small fields of view, that shouldn’t be a problem.

If they are WSI… maybe @petebankhead will have an idea.

Update for M11, also changes label text now to indicate completion of a particular R^2 value.

Still running into errors with this script with certain image types including the LuCa image, while it works for other images with similar bit types. Not sure what is going wrong, but I get a null pointer exception when accessing the pixel values. Works fine for Vectra OMETIFFs and CZIs.

https://gist.github.com/Svidro/68dd668af64ad91b2f76022015dd8a45#file-colocalization-of-channels-per-detection-0-2-0-groovy