class02_ppi_R_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing the Pathway Commons 2 (PC2) database SIF file\n",
    "## CS446/546 class session 2\n",
    "\n",
    "### Goal: count the number of different types of biological interactions in PC2\n",
    "### Approach: retrieve compressed tab-delimited \"edge-list\" file and tabulate \"interaction\" column\n",
    "\n",
    "### Information you will need:\n",
    "\n",
    "- The URL is: http://www.pathwaycommons.org/archives/PC2/v9/PathwayCommons9.All.hgnc.sif.gz\n",
    "- You'll be using the Python modules `gzip`, `timeit`, `pandas`, `urllib.request`, `collections` and `operator`\n",
    "\n",
    "### Other stuff you should do:\n",
    "- Print the first six lines of the uncompressed data file\n",
    "- Use a timer to time how long your program takes\n",
    "- Count how many rows there are in the data file \n",
    "- Estimate the number of proteins in the database; we'll define them operationally as strings in column 1 or column 3, for which the content of column 2 is one of these interactions: 'interacts-with', 'in-complex-with', 'neighbor-of'\n",
    "- Count the total number of unique pairs of interacting molecules (ignoring interaction type)\n",
    "- Count the number rows for each type of interaction in the database\n",
    "- R aficionados:  can you do it using cURL and gunzip, or download.file?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step-by-step instructions for R:\n",
    "\n",
    "- Open a file object representing a stream of the remote, compressed data file, using `urlopen`\n",
    "- Open a file object representing a stream of the uncompressed data file, using `gzip.GzipFile`\n",
    "- Start the timer\n",
    "- Read one line at a time, until the end of the file\n",
    "- Split line on \"\\t\" and pull out the tuple of species1, interaction_type, species2 from the line of text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Class Session 2\n",
    "\n",
    "# this is the URL of the SIF file at Pathway Commons\n",
    "sif_file_url <- \"http://www.pathwaycommons.org/archives/PC2/v9/PathwayCommons9.All.hgnc.sif.gz\"\n",
    "\n",
    "# for starters, we only want three possible interaction types\n",
    "interaction_types_ppi <- c(\"interacts-with\",\"in-complex-with\",\"neighbor-of\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the remote compressed file into a data frame using `readr::read_delim_chunked`; time it using `system.time`; use `head` to print the first six lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Parsed with column specification:\n",
      "cols(\n",
      "  species1 = col_character(),\n",
      "  interaction_type = col_character(),\n",
      "  species2 = col_character()\n",
      ")\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  1.484   0.052   5.462 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>species1</th><th scope=col>interaction_type</th><th scope=col>species2</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>ABCC6         </td></tr>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>ACE2          </td></tr>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>ADAM10        </td></tr>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>ADAM17        </td></tr>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>ADAM9         </td></tr>\n",
       "\t<tr><td>A1BG          </td><td>interacts-with</td><td>AGO1          </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " species1 & interaction\\_type & species2\\\\\n",
       "\\hline\n",
       "\t A1BG           & interacts-with & ABCC6         \\\\\n",
       "\t A1BG           & interacts-with & ACE2          \\\\\n",
       "\t A1BG           & interacts-with & ADAM10        \\\\\n",
       "\t A1BG           & interacts-with & ADAM17        \\\\\n",
       "\t A1BG           & interacts-with & ADAM9         \\\\\n",
       "\t A1BG           & interacts-with & AGO1          \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "species1 | interaction_type | species2 | \n",
       "|---|---|---|---|---|---|\n",
       "| A1BG           | interacts-with | ABCC6          | \n",
       "| A1BG           | interacts-with | ACE2           | \n",
       "| A1BG           | interacts-with | ADAM10         | \n",
       "| A1BG           | interacts-with | ADAM17         | \n",
       "| A1BG           | interacts-with | ADAM9          | \n",
       "| A1BG           | interacts-with | AGO1           | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  species1 interaction_type species2\n",
       "1 A1BG     interacts-with   ABCC6   \n",
       "2 A1BG     interacts-with   ACE2    \n",
       "3 A1BG     interacts-with   ADAM10  \n",
       "4 A1BG     interacts-with   ADAM17  \n",
       "5 A1BG     interacts-with   ADAM9   \n",
       "6 A1BG     interacts-with   AGO1    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# do the filtering while reading the data; reduce disk space and memory usage\n",
    "system.time(interactions_df <- read_delim_chunked(sif_file_url,\n",
    "                                     callback=DataFrameCallback$new(function(df_chunk, pos){ \n",
    "                                        subset(df_chunk, interaction_type %in% interaction_types_ppi)}),\n",
    "                                     chunk_size=10000,\n",
    "                                     delim=\"\\t\",\n",
    "                                     quote=\"\",\n",
    "                                     comment=\"\",\n",
    "                                     col_names=c(\"species1\",\"interaction_type\",\"species2\"),\n",
    "                                     progress=FALSE))\n",
    "\n",
    "# sanity check the resulting data frame\n",
    "head(interactions_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count the number of rows that correspond to protein-protein interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "523498"
      ],
      "text/latex": [
       "523498"
      ],
      "text/markdown": [
       "523498"
      ],
      "text/plain": [
       "[1] 523498"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# how many rows are there in the protein-protein interaction edge-list?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count the number of proteins that participate in protein-protein interactions, using `unique`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "17020"
      ],
      "text/latex": [
       "17020"
      ],
      "text/markdown": [
       "17020"
      ],
      "text/plain": [
       "[1] 17020"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# how many unique proteins are there in the interaction network?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count the number of unique interacting protein-protein pairs, regardless of interaction type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "491784"
      ],
      "text/latex": [
       "491784"
      ],
      "text/markdown": [
       "491784"
      ],
      "text/plain": [
       "[1] 491784"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# how many unique interacting protein pairs are there in the interaction network?\n",
    "length(unique(apply(interactions_df[,c(1,3)], 1, function(my_pair) {\n",
    "    paste(c(min(my_pair), \"-\", max(my_pair)),collapse=\"\")\n",
    "})))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `table` to count the number of each type of interaction in the database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  2.820   0.152   4.315 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                interacts-with                in-complex-with \n",
      "                        369895                         153603 \n",
      "              chemical-affects             catalysis-precedes \n",
      "                        135268                         120948 \n",
      "        controls-expression-of       controls-state-change-of \n",
      "                        110013                         106156 \n",
      "        controls-production-of      consumption-controlled-by \n",
      "                         18482                          16816 \n",
      "   controls-phosphorylation-of                used-to-produce \n",
      "                         15636                          13705 \n",
      "         controls-transport-of                    reacts-with \n",
      "                          6960                           3607 \n",
      "controls-transport-of-chemical \n",
      "                          2847 \n"
     ]
    }
   ],
   "source": [
    "# really quick-and-dirty approach using cURL and gunzip\n",
    "system.time({\n",
    "    system(paste(c(\"curl -s --compressed \", sif_file_url, \" | gunzip > pc.sif\"), collapse=\"\"), \n",
    "       intern=TRUE)\n",
    "my_df <- read.table(file=\"pc.sif\",\n",
    "                    sep=\"\\t\",\n",
    "                    comment.char=\"\",\n",
    "                    quote=\"\",\n",
    "                    col.names=c(\"species1\",\n",
    "                                \"interaction_type\",\n",
    "                                \"species2\"),\n",
    "                    stringsAsFactors=FALSE)\n",
    "interactions_df5 <- subset(my_df,\n",
    "                          interaction_type %in% interaction_types_ppi)})\n",
    "\n",
    "all(interactions_df5 == interactions_df)\n",
    "\n",
    "\n",
    "## TODO:  make a factor out of the interaction_type column, call table on it, and then sort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       " 26.148   0.212  26.795 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# alternative approach which is kind of slow and (transiently) memory-hungry; read the entire file into memory and \n",
    "# then process the text contents of the file, line by line\n",
    "\n",
    "system.time({\n",
    "alltext <- readLines(gzcon(url(sif_file_url)))\n",
    "nlines <- length(alltext)\n",
    "interactions_df4 <- data.frame(do.call(rbind, lapply(1:length(alltext),\n",
    "                                         function(i) {\n",
    "                                             split_line <- strsplit(alltext[i], \"\\t\")[[1]]\n",
    "                                             names(split_line) <- c(\"species1\",\"interaction_type\",\"species2\")\n",
    "                                             split_line\n",
    "                                         })))\n",
    "interactions_df4 <- subset(interactions_df4, interaction_type %in% interaction_types_ppi)\n",
    "})\n",
    "            \n",
    "# sanity check that this slow approach at least gave us correct results\n",
    "all(interactions_df4 == interactions_df)\n",
    "rm(interactions_df4)\n",
    "rm(alltext)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  2.988   0.000   4.873 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# quick and dirty approach which doesn't require looking up any newfangled commands, but transiently eats up\n",
    "# disk space and memory\n",
    "fname <- tempfile()\n",
    "gzfname <- paste(c(fname, \".gz\"),collapse=\"\")\n",
    "system.time({\n",
    "   download.file(sif_file_url, destfile=gzfname)\n",
    "   interactions_df3 <- subset(read.table(file=gzfile(gzfname),\n",
    "                              sep=\"\\t\",\n",
    "                              comment.char=\"\",\n",
    "                              quote=\"\",\n",
    "                              col.names=c(\"species1\",\"interaction_type\", \"species2\"),\n",
    "                              stringsAsFactors=FALSE),\n",
    "                           interaction_type %in% interaction_types_ppi)})\n",
    "unlink(gzfname)\n",
    "\n",
    "# sanity check the results\n",
    "all(interactions_df3 == interactions_df)\n",
    "\n",
    "# clean up\n",
    "rm(interactions_df3)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}