Workshop Workflow PAX6-unsolved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# InterMineR Workshop Use Case\n",
    "\n",
    "\n",
    "\n",
    "We are going to re-create the workflow we did using the web interface using the R API.\n",
    "\n",
    "The basic steps are:  \n",
    "\n",
    "1. Load the InterMine library and choose an InterMine to query.\n",
    "2. Query 1: Diabetes Genes: Fetch a list of genes that are associated with diabetes\n",
    "3. Query 2: PAX6 + Pancreas: Fetch a list of genes that have medium or high expression in the pancreas and are in our PAX6 targets list\n",
    "4. Intersection: Find which genes are present in _both_ Query 1 and Query2.\n",
    "5. GWAS: Compare the intersection of the previous query with results from GWAS studies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting started - Load InterMineR and choose an InterMine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the InterMine library. If it's not already installed, visit\n",
    "https://bioconductor.org/packages/release/bioc/html/InterMineR.html and follow the instructions to install.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(InterMineR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to query human data - so let's look and see what InterMines are available.  \n",
    "The method we want to use is `listMines()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "listMines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, let's select HumanMine from the list, and store it in a variable called `humanMine`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, now let's tell InterMineR that we want to use HumanMine for our queries. Store it in a variable called `im`.\n",
    "\n",
    "**Important:** you'll need an API token for this part so you can access your HumanMine account. You can get your token by logging into [HumanMine](http://www.humanmine.org/) and going to the account details tab within MyMine. Cut and paste your token into the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the syntax to initialise an InterMine is as follows:\n",
    "# im <- initInterMine(mine=SOME-INTERMINE, \"YOUR TOKEN HERE\")\n",
    "im <- "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Query: Diabetes Genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our first query will be to select all human genes that are associate with diabetes.   \n",
    "\n",
    "This will require two **constraints**:\n",
    "1. Ensure all genes are human genes. Make sure `Gene.organism.name` is equal (use the `=` operator) to `Homo sapiens`. (HumanMine contains some non-human genes for homology query purposes)\n",
    "2. Restrict results to genes that are associated with `diabetes`. This constraint will be asking for  `Gene.diseases.name` to contain (use the operator `CONTAINS`) the string `\"diabetes\"`\n",
    "\n",
    "We'll also **select two columns to view**:  \n",
    "1. `Gene.primaryIdentifier`\n",
    "2. `Gene.symbol`\n",
    "\n",
    "Store the first query in a variable valled `query1Diabetes`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll be using the setQuery method here, with two arguments: select and where.\n",
    "# \n",
    "# Syntax: \n",
    "# setQuery( \n",
    "#  select = c(\"VIEW.NAME1\", \"VIEW.NAME2\", .... ),\n",
    "#  where = setConstraints(\n",
    "#    paths = c(\"CONSTRAINT.PATH.A\", \"CONSTRAINT.PATH.B\"),\n",
    "#    operators = c(\"CONSTRAINT.OPERATOR.A\", \"CONSTRAINT.OPERATOR.B\"),\n",
    "#    values = list(\"CONSTRAINT.VALUE.A\",\"CONSTRAINT.VALUE.B\")\n",
    "#  )\n",
    "#)\n",
    "\n",
    "\n",
    "query1Diabetes <- setQuery( \n",
    "  # here we're choosing which columns of data we'd like to see\n",
    "  select = c(  ),\n",
    "    # set the logic for constraints. The first constraint is the first path+operator+value, \n",
    "    # e.g. Gene.organism.name = Homo sapiens, and the second constraint is the combination \n",
    "    # of the second path+operator+value, e.g. Gene.diseases.name CONTAINS diabetes\n",
    "  where = setConstraints(\n",
    "    paths = c(  ),\n",
    "    operators = c(  ),\n",
    "    values = list(  )\n",
    "  )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question to ponder:** why did we use `=` for our Homo sapiens constraint, but `CONTAINS` for our diabetes constraint?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anyway, we've set the query up, so now let's actually run it. Use the method `runQuery` and store the results in a variable called `query1DiabetesResults`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Syntax: runQuery(im,QUERY_TO_RUN)\n",
    "\n",
    "query1DiabetesResults <- #put the code to run the query here\n",
    "\n",
    "# and let's print out the first few results to make sure it looks like we'd expect:\n",
    "# Syntax: head(THING_TO_PRINT) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query 2: Pax6 targets that have high expression in the Pancreas\n",
    "\n",
    "This time we're creating another query, but with slightly more complex constraints. We're looking for genes that are in the public HumanMine list `PL_Pax6_Targets`, that are also expressed in the pancreas at a `High` or `Medium` level. \n",
    "\n",
    "We'll need a few more **constraints** than we did in Query 1:  \n",
    "1. all `Gene`s should be `IN` the list `PL_Pax6_Targets`\n",
    "2. `Gene.proteinAtlasExpression.level` should be set to `High` OR `Medium`. This will require two constraints, one for each of medium and high.\n",
    "3. `Gene.proteinAtlasExpression.tissue.name` should be equal to `Pancreas`\n",
    "\n",
    "We'd also like to see a few more **columns** this time: \n",
    "1. The `Gene`'s `primaryIdentifier` and `symbol`\n",
    "2. The following expression data from Protein Atlas: \n",
    "    - `Gene.proteinAtlasExpression.cellType`\n",
    "    - `Gene.proteinAtlasExpression.level`\n",
    "    - `Gene.proteinAtlasExpression.tissue.name`\n",
    "    \n",
    "Store it in a variable called `query2UpInPancreas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We don't want to see *all* genes and their expression. \n",
    "# Let's narrow it down a little by constraining it to genes that are \n",
    "# of interest, using setConstraints. \n",
    "# Syntax: \n",
    "# \n",
    "# setConstraints(\n",
    "#    paths = c(\"CONSTRAINT.PATH.A\", \"CONSTRAINT.PATH.B\"),\n",
    "#    operators = c(\"CONSTRAINT.OPERATOR.A\", \"CONSTRAINT.OPERATOR.B\"),\n",
    "#    values = list(\"CONSTRAINT.VALUE.A\",\"CONSTRAINT.VALUE.B\")\n",
    "# )\n",
    "\n",
    "query2UpInPancreasConstraint <- setConstraints(\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excellent - we've defined the constraints we want. We still need to choose which columns to view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new query.\n",
    "# Syntax: \n",
    "# newQuery(\n",
    "#  view = c(\"VIEW.NAME.1\", \"VIEW.NAME.2\"), \n",
    "#  constraintLogic = \"X and (Y or Z)\"      #constraintLogic is optional.\n",
    "# )    \n",
    "query2UpInPancreas = newQuery(\n",
    "  # Choose which columns of data we'd like to see\n",
    "  view = c(\n",
    "      \n",
    "  ),\n",
    "  # set the logic for constraints. This means our pancreas expression level \n",
    "  # is EITHER Medium (B) or High (C), but not both.\n",
    "  # --\n",
    "  # Note: Constraint logic only needs to be set if you wish to use OR. All other\n",
    "  # constraints have AND logic applied by default. \n",
    "  constraintLogic = # Add your logic here\n",
    ")\n",
    "\n",
    "# Add the constraint to our expressed pancreas query (previously we just _defined_ the constraint)\n",
    "# Syntax: QUERYNAME$where <- CONSTRAINT\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, that was just setting up the query - we haven't run it yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we have the query set up the way we want, let's actually *run* the query! \n",
    "# Syntax: runQuery(im,QUERY_TO_RUN)\n",
    "query2UpInPancreasResults \n",
    "\n",
    "# Show me the first few results please! \n",
    "head(query2UpInPancreasResults) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intersection: Which genes overlap in Query1 and Query2?\n",
    "\n",
    "Let's check which genes are in BOTH lists that we've created. To do this, we'll strip down the columns we have to just primary identifiers, and then run a list intersect function.\n",
    "\n",
    "1. Extract the `\"Gene.primaryIdentifier\"` from both of your result lists - `query1DiabetesResults` and `query2UpInPancreasResults`.\n",
    "2. Run `intersect(col1, col2)` on the columns you extracted\n",
    "3. Print the results and admire them! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the primaryIdentifier columns from query1 (diabetes genes) and query 2 (upexpressed in pancreas)\n",
    "# \n",
    "# Syntax to get an individual column view: \n",
    "# data.frame[[\"COLUMN.NAME\"]]\n",
    "primaryIdentifiersDiabetes <- # put your code here\n",
    "primaryIdentifiersPancreas <- # put your code here\n",
    "\n",
    "# Find the intersection of the two vectors of primary identifiers, and store it in a variable\n",
    "#\n",
    "# Syntax: intersect(VECTOR_1, VECTOR_2)\n",
    "\n",
    "# Show the results using the print(VARIABLE_TO_PRINT) method. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GWAS: Compare the intersection above with results from GWAS studies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we fed the intersected list from above back into another query to see if there was any association of these genes with diabetes phenotypes according to GWAS studies. Note that we now start our query from the `GWAS` class.\n",
    "\n",
    "Start by making a constraint with the following values:\n",
    "1. `\"GWAS.results.pValue\"` is less than or equals (use the operator `<=`) to `\"1e-04\"`. (1e-04 is the same as 0.0001).\n",
    "2. `\"GWAS.results.phenotype\"` must contain (operator `CONTAINS`) the word `diabetes`.\n",
    "3. `\"GWAS.results.associatedGenes.primaryIdentifier\"` should be equal to each of the values from your previous step. You'll need to repeat this constraint once per identifier. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we set up the constraints. The last three constraints are the \n",
    "# diabetesAndPancreas result genes from our last query.\n",
    "#\n",
    "# Syntax: \n",
    "# \n",
    "# setConstraints(\n",
    "#    paths = c(\"CONSTRAINT.PATH.A\", \"CONSTRAINT.PATH.B\"),\n",
    "#    operators = c(\"CONSTRAINT.OPERATOR.A\", \"CONSTRAINT.OPERATOR.B\"),\n",
    "#    values = list(\"CONSTRAINT.VALUE.A\",\"CONSTRAINT.VALUE.B\")\n",
    "# )\n",
    "query3GWASConstraints <- setConstraints(\n",
    "    # put your constraints here\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we've set our constraints up nicely, let's choose which columns we want to view. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query3GWAS <- newQuery( \n",
    "  # Quite a few columns this time!\n",
    "  view = c(\"GWAS.results.associatedGenes.primaryIdentifier\",\n",
    "    \"GWAS.results.associatedGenes.symbol\", \"GWAS.results.associatedGenes.name\",\n",
    "    \"GWAS.results.SNP.primaryIdentifier\", \"GWAS.results.pValue\", \"GWAS.results.phenotype\",\n",
    "    \"GWAS.firstAuthor\", \"GWAS.name\", \"GWAS.publication.pubMedId\",\n",
    "    \"GWAS.results.associatedGenes.organism.shortName\"),\n",
    "    # set the logic for constraints. Remember that we want our results\n",
    "    # to include any one of the three genes we found in the list of diabetes+pancreas genes\n",
    "    # so we need to use some OR logic.\n",
    "  constraintLogic = #Fill the constraint logic in here. Remember the syntax: \"A and (B or C)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the constraints to the query, and then run it - we've filled this step in for you"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add constraint\n",
    "query3GWAS$where <- query3GWASConstraints\n",
    "#run query\n",
    "query3GWASResults <- runQuery(im, query3GWAS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's view those results..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query3GWASResults"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's take a look at the unique gene symbols that were returned:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GWASIds <- query3GWASResults[\"GWAS.results.associatedGenes.symbol\"]\n",
    "unique(GWASIds)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}