IDR0021-Segmentation.Rmd

## Image Segmentation using the idr0021 dataset

This notebook uses the idr0021 dataset ([idr0021-lawo-pericentriolarmaterial/experimentA](https://idr.openmicroscopy.org/webclient/?show=project-51)) and tries to reproduce some of the analysis published in ['Subdiffraction imaging of centrosomes reveals higher-order organizational features of pericentriolar material'](https://doi.org/10.1038/ncb2591); in particular to create a figure similar to [Figure 1](https://www.nature.com/articles/ncb2591/figures/1) of the article.

### Load libraries

```{r}
library(romero.gateway)
library(EBImage, warn.conflicts = FALSE)
```

### Connect to server

```{r}
server <- OMEROServer(host = 'wss://idr.openmicroscopy.org/omero-ws', port = 443L, username='public', password='public')
server <- connect(server)
paste('Successfully logged in as', server@user$getUserName())
```

### Manual segmentation

```{r}
# Load Image 'siControl_N20_Cep215_I_20110411_Mon-1509_0_SIR_PRJ.dv' from dataset 'CDK5RAP2-C'
imageId <- 1884839 
image <- loadObject(server, "ImageData", imageId)
paste("Image", imageId, "loaded.")
```

#### Load pixel values and display image

```{r}
# There is just one plane, so z = 1 and t = 1
z <- 1
t <- 1

# Load the second channel
channelIndex <- 2

pixels <- getPixelValues(image, z, t, channelIndex)

ebimage <- EBImage::Image(data = pixels, colormode = 'Grayscale')
ebimage <- normalize(ebimage)
EBImage::display(ebimage)
```

#### Threshold

Feel free to play around with the parameters to get a better result.

```{r}
# Threshold
threshImg <- thresh(ebimage, w=15, h=15, offset=0.1)

# Remove noise
threshImg <- medianFilter(threshImg, size=3)

# Fill holes
threshImg <- fillHull(threshImg)          

# Distinguish objects
threshImg <- bwlabel(threshImg)

# Show the segmented image
EBImage::display(colorLabels(threshImg))
```

#### Compute features (measurements)

```{r}
shapes = computeFeatures.shape(threshImg)
moments = computeFeatures.moment(threshImg)

# merge the two dataframes together into one 'features' dataframe
features <- merge(shapes, moments, by=0, all=TRUE)
features
```

#### Save ROIs to OMERO

If you'd be connected to an OMERO server on which you have write permissions, you could save the segmentation result as ROIs on the image. This will not work against the public IDR OMERO server. However, that's how you would do it:

```{r}
#' Create ROIs from a features table
#'
#' @param features The shape and moments generated by EBImage computeFeatures
#' @return A dataframe specifying the x, y, rx, ry, w, h, theta and text parameters of the ROIs
createROIs <- function(features) {
    rois <- data.frame(x = c(0), y = c(0), rx = c(0), ry = c(0), w = c(0), h = c(0), theta = c(0), text = c('remove'), stringsAsFactors = FALSE)
    for (index in 1:length(features[,1])) {
        x <- features[index,8]
        y <- features[index,9]
        r1 <- features[index,10]
        ecc <- features[index,11]
        r2 <- sqrt(r1^2 * (1 - ecc^2))
        theta <- features[index,12]
        rois <- rbind(rois, c(x, y, r1, r2, NA, NA, -theta, as.character(index)))
    }
    rois <- rois[-1,]
    rownames(rois) <- c()
    return(rois)
}

rois <- createROIs(features)
rois

# addROIs creates an ROI for each entry in the dataframe specified by the 'coords' parameter
# addROIs(image, coords = rois) # can't run this against public IDR server

# Save the features dataframe on the image:
# attachDataframe(image, features, "ROI features")) # again this needs write permission

# Alternatively add the features dataframe as csv file to the image:
# csvFile <- "/tmp/ROI_features.csv" 
# write.csv(features, file = csvFile)
# fileannotation <- attachFile(image, csvFile) # again this needs write permission
```

### Automate analysis

```{r}
#' Performs an image segmentation to find the largest ROI of the image
#' and returns some features of interest (area, perimeter and diameter).
#' Optionally: Creates an ROI for each object detected by the segmentation.
#'
#' @param image The image to segment
#' @param channelName The channel to be taken into account
#' @param df The dataframe to add the features to (channelName, imageName, ImageID, ROIIndex, area, perimeter, diameter)
#' @param saveROIs Flag if ROIs should be created and attached to the image (default: FALSE)
#' @return The dataframe with the features
analyzeImage <- function(image, channelName, df, saveROIs = FALSE) {
    # Find the channel index
    chnames <- getChannelNames(image)
    chindex <- match(channelName, chnames, nomatch = 0)
    if (chindex == 0) {
      message (paste("Could not resolve channel name, skipping ", image@dataobject$getId()))
      return(df)
    }
    
    # Load the pixels
    pixels <- getPixelValues(image, 1, 1, chindex)
    ebi <- EBImage::Image(data = pixels, colormode = 'Grayscale')
    ebi <- normalize(ebi)
     
    # this is our segmentation workflow from above
    threshImg <- thresh(ebi, w=15, h=15, offset=0.1)
    threshImg <- medianFilter(threshImg, size=3)
    threshImg <- fillHull(threshImg)          
    threshImg <- bwlabel(threshImg)
    
    # Calculate the features
    shapes = suppressMessages(computeFeatures.shape(threshImg))
    moments = suppressMessages(computeFeatures.moment(threshImg))
    features <- merge(shapes, moments, by=0, all=TRUE)
    
    if (length(features[,1])>1) {
        # Add the ROIs to the image
        if (saveROIs) {
            rois <- createROIs(features)
            addROIs(image, coords = rois)
        }
        
        # Add the interesting properties (area, perimeter and diameter)
        # of the largest object together with channel name, image name, image id 
        # and roi index to the dataframe
        features <- features[order(-features[,2]),]
        diameter <- features[1,4]*2
        df <- rbind(df, c(channelName, image@dataobject$getName(), image@dataobject$getId(), features[1,1], features[1,2], features[1,3], diameter))
    }
    return(df)
}
```

#### Run analysis on the a whole dataset

This will analyse all images of the 'CDK5RAP2-C' dataset.

```{r}
datasetId <- 51

channelName <- 'CDK5RAP2-C'

dataset <- loadObject(server, "DatasetData", datasetId)

# Keep the channel name, image name, image id, area, perimeter, and diameter of the largest ROIs
result <- data.frame(Channel = c('remove'), ImageName = c('remove'), Image = c(0), ROIIndex = c(0), Area = c(0), Perimeter = c(0), Diameter = c(0), stringsAsFactors = FALSE)

images <- getImages(dataset)
for (image in images) {
    result <- tryCatch({
        analyzeImage(image, channelName, result, saveROIs = TRUE)
    }, warning = function(war) {
        message(paste("WARNING:  ", image@dataobject$getId(),war))
        return(result)
    }, error = function(err) {
        message(paste("ERROR:  ", image@dataobject$getId() ,err))
        return(result)
    }, finally = {
    })
}

result <- result[-1,]
rownames(result) <- c()

# set the correct datatypes
result$Channel <- as.factor(result$Channel)
result$Area <- as.numeric(result$Area)
result$Perimeter <- as.numeric(result$Perimeter)
result$Diameter <- as.numeric(result$Diameter)

# set the correct pixel size
pxSizeInNM <- 40
result$Area <- result$Area * pxSizeInNM ^ 2
result$Perimeter <- result$Perimeter * pxSizeInNM
result$Diameter <- result$Diameter * pxSizeInNM

result
```

### Statistical analysis of the whole idr0021 dataset

Running the segmentation over the whole project takes some time. This has already been done, lets take the short cut and load the data:

```{r}
load("IDR0021.RData") # load 'centrioles' dataframe
centrioles
```

#### Plot the data

```{r}
# Order the datasets ascending by mean diameter 
ag <-aggregate(centrioles$Diameter ~ centrioles$Dataset, centrioles, median)
ordered <- factor(centrioles$Dataset, levels=ag[order(ag$`centrioles$Diameter`), 'centrioles$Dataset'])

# Draw the plot 
plot(centrioles$Diameter ~ ordered, ylab='Diameter (nm)', xlab="Protein", cex.axis=0.5)
```

#### One way ANOVA

```{r}
fit <- aov(centrioles$Diameter ~ centrioles$Dataset)
summary(fit)
```

#### Two-sample Wilcoxon test of all pairwise combinations

```{r}
# Two-sample Wilcoxon test ('Mann-Whitney') of all pairwise combinations:
combins <- combn(levels(centrioles$Dataset), 2)
params_list <- split(as.vector(combins), rep(1:ncol(combins), each = nrow(combins)))
testResults <- data.frame()
for (param in params_list) {
  testdf <- subset(centrioles, centrioles$Dataset %in% param)
  pval <- wilcox.test(formula = Diameter ~ Dataset, data = testdf)$p.value
  testResults <- rbind(testResults, data.frame(Protein_1=param[1], Protein_2=param[2], p_value=pval))
}
testResults
```