R/idr0021_rois.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IDR0021 - ROIs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple example which reads the 'batch_roi_export' table, performs some basic statistics on it and creates a plot (which should be similar to [Figure 1](https://www.nature.com/articles/ncb2591/figures/1) in the paper ['Subdiffraction imaging of centrosomes reveals higher-order organizational features of pericentriolar material'](https://doi.org/10.1038/ncb2591)). The plot and results are then uploaded to the server and attached to the project."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the romero.gateway library and connect to the server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(warn=-1) # Disable warnings to suppress unnecessary output\n",
    "\n",
    "library(romero.gateway)\n",
    "\n",
    "user_name = readline('Username: ')\n",
    "user_password <- getPass::getPass('OMERO password: ')\n",
    "server <- OMEROServer(host = 'workshop.openmicroscopy.org', port = 4064L, username = user_name, password = user_password)\n",
    "server <- connect(server)\n",
    "paste('Successfully logged in as', server@user$getUserName())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data from the 'batch_roi_export' table as R dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the project and table file\n",
    "projectName <- 'idr0021'\n",
    "tableName <- 'Summary_from_Fiji'\n",
    "\n",
    "searchResult <- searchFor(server, Project, 'Name', projectName)\n",
    "project <- searchResult[[1]]\n",
    "\n",
    "annos <- getAnnotations(server, 'ProjectData', getOMEROID(project), nameFilter = tableName)\n",
    "annotationFileID = as.integer(annos$FileID)\n",
    "\n",
    "\n",
    "# Load the table (column 14 and 15, Dataset and bouding_box) directly as R dataframe\n",
    "df <- loadDataframe(project, annotationFileID, columns=c(14, 15))\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CEP120 is represented by two datasets'CEP120/20111106' and 'CEP120/20111209',\n",
    "# just combine them to 'CEP120'\n",
    "df$Dataset <- as.character(df$Dataset)\n",
    "df[ df == 'CEP120/20111106' ] <- 'CEP120'\n",
    "df[ df == 'CEP120/20111209' ] <- 'CEP120'\n",
    "df$Dataset <- as.factor(df$Dataset)\n",
    "\n",
    "\n",
    "# Use the bounding box of the biggest shapes for the specified channel\n",
    "df$diameter <- sqrt(df$Bounding_Box)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Order the data from lowest to highest mean diameter (to match the order of figure 1)\n",
    "ag <-aggregate(df$diameter ~ df$Dataset, df, mean)\n",
    "orderedDatasets <- factor(df$Dataset, levels=ag[order(ag$`df$diameter`), 'df$Dataset'])\n",
    "\n",
    "plot(df$diameter ~ orderedDatasets, ylab='Bounding Box', xlab=\"Protein\", cex.axis=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Is there a significant difference between the proteins?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# One-way analysis of variance:\n",
    "fit <- aov(df$diameter ~ df$Dataset)\n",
    "summary(fit)\n",
    "\n",
    "# Two-sample Wilcoxon test ('Mann-Whitney') of all pairwise combinations:\n",
    "combins <- combn(levels(df$Dataset), 2)\n",
    "params_list <- split(as.vector(combins), rep(1:ncol(combins), each = nrow(combins)))\n",
    "testResults <- data.frame()\n",
    "for (param in params_list) {\n",
    "  testdf <- subset(df, df$Dataset %in% param)\n",
    "  pval <- wilcox.test(formula = diameter ~ Dataset, data = testdf)$p.value\n",
    "  testResults<-rbind(testResults, data.frame(Protein_1=param[1], Protein_2=param[2], p_value=pval))\n",
    "}\n",
    "testResults"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Attach the results to OMERO project and close connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The plot:\n",
    "tmpfile <- \"/tmp/Bounding_Box_By_Protein.png\"\n",
    "png(tmpfile, width = 4, height= 2.5, units = \"in\", res = 300, pointsize = 6)\n",
    "plot(df$diameter ~ orderedDatasets, ylab='Bounding Box', xlab=\"Protein\", cex.axis=0.5)\n",
    "dev.off()\n",
    "invisible(attachFile(project, tmpfile)) # Wrapped in invisible(...) only to suppress unnecessary output\n",
    "\n",
    "\n",
    "# The results from the Mann-Whitney test as OMERO.table (HDF)\n",
    "invisible(attachDataframe(project, testResults, \"Mann-Whitney-pValues\"))\n",
    "\n",
    "# and as CSV file:\n",
    "tmpfile <- \"/tmp/Mann-Whitney-pValues.csv\"\n",
    "write.csv(testResults, file = tmpfile)\n",
    "invisible(attachFile(project, tmpfile)) \n",
    "\n",
    "# Finally disconnect again\n",
    "invisible(disconnect(server))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}