class01_ppi_R_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing the Pathway Commons 2 (PC2) database SIF file\n",
    "## CS446/546 class session 1\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",
    "### What is the \".sif\" file format?\n",
    "\n",
    "SIF stands for Simple Interaction File. The format is like this:\n",
    "```\n",
    "A1BG    controls-expression-of  A2M\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",
    "...\n",
    "```\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": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Parsed with column specification:\n",
      "cols(\n",
      "  species1 = \u001b[31mcol_character()\u001b[39m,\n",
      "  interaction_type = \u001b[31mcol_character()\u001b[39m,\n",
      "  species2 = \u001b[31mcol_character()\u001b[39m\n",
      ")\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  3.013   0.182   4.179 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A tibble: 6 × 3</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>species1</th><th scope=col>interaction_type</th><th scope=col>species2</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
       "</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": [
       "A tibble: 6 × 3\n",
       "\\begin{tabular}{r|lll}\n",
       " species1 & interaction\\_type & species2\\\\\n",
       " <chr> & <chr> & <chr>\\\\\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",
       "A tibble: 6 × 3\n",
       "\n",
       "| species1 &lt;chr&gt; | interaction_type &lt;chr&gt; | species2 &lt;chr&gt; |\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"
      ],
      "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": [
    "library(readr)\n",
    "# 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": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "508480"
      ],
      "text/latex": [
       "508480"
      ],
      "text/markdown": [
       "508480"
      ],
      "text/plain": [
       "[1] 508480"
      ]
     },
     "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "17531"
      ],
      "text/latex": [
       "17531"
      ],
      "text/markdown": [
       "17531"
      ],
      "text/plain": [
       "[1] 17531"
      ]
     },
     "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;\n",
    "use `apply` to generate a unique key for each pair of proteins, then use the `unique` function on the keys and then the `length` function to compute the number of unique keys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "475553"
      ],
      "text/latex": [
       "475553"
      ],
      "text/markdown": [
       "475553"
      ],
      "text/plain": [
       "[1] 475553"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# how many unique interacting protein pairs are there in the \n",
    "# interaction network?\n",
    "\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": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  2.287   0.109   4.778 "
      ]
     },
     "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",
      "              chemical-affects                 interacts-with \n",
      "                        492765                         325616 \n",
      "               in-complex-with       controls-state-change-of \n",
      "                        182864                         182450 \n",
      "            catalysis-precedes         controls-expression-of \n",
      "                        149013                         123232 \n",
      "     consumption-controlled-by         controls-production-of \n",
      "                         22830                          21494 \n",
      "   controls-phosphorylation-of                used-to-produce \n",
      "                         17029                          14486 \n",
      "         controls-transport-of                    reacts-with \n",
      "                          7574                           3927 \n",
      "controls-transport-of-chemical \n",
      "                          3322 \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",
    "print(sort(table(as.factor(my_df$interaction_type)), decreasing=TRUE))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       " 26.335   0.287  27.388 "
      ]
     },
     "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": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  1.808   0.074   4.484 "
      ]
     },
     "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 \n",
    "# newfangled commands, but transiently eats up 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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}