This notebook uses the idr0021 dataset (idr0021-lawo-pericentriolarmaterial/experimentA) and tries to reproduce some of the analysis published in 'Subdiffraction imaging of centrosomes reveals higher-order organizational features of pericentriolar material'; in particular to create a figure similar to Figure 1 of the article.
library(romero.gateway) library(EBImage, warn.conflicts = FALSE)
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())
# 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.")
# 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)
Feel free to play around with the parameters to get a better result.
# 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))
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
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:
#' 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
#' 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) }
This will analyse all images of the 'CDK5RAP2-C' dataset.
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
Running the segmentation over the whole project takes some time. This has already been done, lets take the short cut and load the data:
load("IDR0021.RData") # load 'centrioles' dataframe centrioles
# 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)
fit <- aov(centrioles$Diameter ~ centrioles$Dataset) summary(fit)
# 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.