EXTRA Challenge exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: right;\" src=\"files/thehyve_logo.png\">\n",
    "# EXTRA Challenge exercise\n",
    "---------------\n",
    "Copyright (c) 2015 The Hyve B.V. This notebook is licensed under the GNU General Public License, version 3. Authors: Ruslan Forostianov, Pieter Lukasse. Parts of the R script: Copyright 2014 Janssen Research & Development, LLC, & Copyright (c) 2015 The Hyve B.V., based on demoTransmartRClientCommands.R also made available as GPL v3 in this Jupyter home dir."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare: authenticate with tranSMART \n",
    "------------------\n",
    "\n",
    "Authenticate with tranSMART first if you want to execute any of the analysis in the boxes below again. \n",
    "\n",
    "**Step 1**: Please open URL http://localhost:8080/transmart/oauth/authorize?response_type=code&client_id=api-client&client_secret=api-client&redirect_uri=http%3A%2F%2Flocalhost%3A8080%2Ftransmart%2Foauth%2Fverify\n",
    "\n",
    "**Step 2**: paste token in the token parameter below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "require(\"transmartRClient\")\n",
    "connectToTransmart(\"http://localhost:8080/transmart\", prefetched.request.token = \"OjNGZi\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the output above is: Authentication completed. TRUE , \n",
    "then you can continue below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1 - Retrieving metadata and clinical data\n",
    "-----------------\n",
    "\n",
    "## Get studies and observations data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "study <- \"GSE8581\"  \n",
    "\n",
    "# Retrieve Clinical Data\n",
    "allObservations <- getObservations(study, as.data.frame = T)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2 - Retrieving molecular (aka \"high dimensional\") data\n",
    "------------------\n",
    "\n",
    "\n",
    "## Downloading the expression data for our chosen study\n",
    "\n",
    "This can take a while (~1 minute)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dataDownloaded <- getHighdimData(study.name = study, concept.match = \"Lung\", projection = \"log_intensity\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "summary(dataDownloaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare the data for easy usage in different standard R functions\n",
    "The steps below show how the table above is processed into a simple table that contains only expression values + an extra feature of having patient identifiers as row names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# select gene expression data, which is the data *excluding* columns 1 to 6:\n",
    "data<-dataDownloaded[[\"data\"]]\n",
    "expression_data<-data[,-c(1:6)]\n",
    "expression_data[1:3,1:3]\n",
    "dim(expression_data)\n",
    "# add patientId as the row name for the expression_data matrix:\n",
    "rownames(expression_data)<-data$patientId\n",
    "expression_data[1:3,1:3]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating the PCA plot \n",
    "If the dimensions of the expression_data table are large, you may want to create a subset of the data first. Here we use a probelist as a *subset* for the probes, based on the list found in: \"Bhattacharya S., Srisuma S., Demeo D. L., et al.,  Molecular biomarkers for quantitative and discrete COPD phenotypes.American Journal of Respiratory Cell and Molecular Biology. 2009;40(3):359–367. doi: 10.1165/rcmb.2008-0114OC.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Select a subset of the probes:\n",
    "probeNames<- c(\"1552622_s_at\",\"1555318_at\",\"1557293_at\",\"1558280_s_at\",\"1558411_at\",\"1558515_at\",\"1559964_at\",\"204284_at\",\"205051_s_at\",\"205528_s_at\",\"208835_s_at\",\"209377_s_at\",\"209815_at\",\"211548_s_at\",\"212179_at\",\"212263_at\",\"213156_at\",\"213269_at\",\"213650_at\",\"213878_at\",\"215359_x_at\",\"215933_s_at\",\"218352_at\",\"218490_s_at\",\"220094_s_at\",\"220906_at\",\"220925_at\",\"222108_at\",\"224711_at\",\"225318_at\",\"225595_at\",\"225835_at\",\"225892_at\",\"226316_at\",\"226492_at\",\"226534_at\",\"226666_at\",\"226800_at\",\"227095_at\",\"227105_at\",\"227148_at\",\"227812_at\",\"227852_at\",\"227930_at\",\"227947_at\",\"228157_at\",\"228630_at\",\"228665_at\",\"228760_at\",\"228850_s_at\",\"228875_at\",\"228963_at\",\"229111_at\",\"229572_at\",\"230142_s_at\",\"230986_at\",\"232014_at\",\"235423_at\",\"235810_at\",\"238712_at\",\"238992_at\",\"239842_x_at\",\"239847_at\",\"241936_x_at\",\"242389_at\")\n",
    "#note: this is because R automatically prepends \"X\" in front of column names that start with a numerical value. Therefore prepend \"X\"\n",
    "probeNames<- paste(\"X\", probeNames, sep = \"\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# select only the cases and controls (excluding the patients for which the lung disease is not specified). Note: in the observation table the database IDs \n",
    "# are used to identify the patients and not the patient IDs that are used in the gene expression dataset\n",
    "cases <- allObservations$observations$subject.id[allObservations$observations$'Subjects_Lung Disease' == \"chronic obstructive pulmonary disease\"]\n",
    "controls <- allObservations$observations$subject.id[allObservations$observations$'Subjects_Lung Disease' == \"control\"]\n",
    "par(pin=c(5,2))\n",
    "barplot(c(length(cases),length(controls)), main=\"Cases and controls\", horiz=TRUE,\n",
    "  names.arg=c(\"cases\", \"controls\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the data\n",
    "\n",
    "Some basic data preps to separate cases and controls and make the data suitable for passing on to the **prcomp()** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# now we have the *internal database* IDs for the patients, but we need to get the patient IDs because \n",
    "# this is the index of the expression_data matrix. \n",
    "# These can be retrieved from the subjectInfo table: \n",
    "subjectInfo <- allObservations$subjectInfo\n",
    "patientIDsCase    <- subjectInfo$subject.inTrialId[ subjectInfo$subject.id %in% cases ] \n",
    "patientIDsControl <- subjectInfo$subject.inTrialId[ subjectInfo$subject.id %in% controls] \n",
    "\n",
    "# patient sets containing case and control patientIDs\n",
    "patientSets <- c(patientIDsCase, patientIDsControl)\n",
    "patientSets <- patientSets[which(patientSets %in% rownames(expression_data))]\n",
    "# make a subset of the data based on the selected patientSets and the probelist, and transpose the \n",
    "# table so that the rows now contain probe names\n",
    "subset<-t(expression_data[patientSets,probeNames]) \n",
    "# for ease of recognition: append \"Case\" and \"Control\" to the patient names\n",
    "colnames(subset)[colnames(subset)%in% patientIDsCase] <- paste(colnames(subset)[colnames(subset)%in% patientIDsCase],\"Case\", sep=\"_\" )\n",
    "colnames(subset)[colnames(subset)%in% patientIDsControl] <- paste( colnames(subset)[colnames(subset)%in% patientIDsControl] , \"Control\",sep= \"_\")\n",
    "\n",
    "# there is one patient that seems to be an outlier: GSE8581GSM212810_Case.\n",
    "# remove this outlier and plot heatmap again\n",
    "subset_without_outlier <- subset[,colnames(subset)!= \"GSE8581GSM212810_Case\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# PCA analysis : \n",
    "options(warn=-1)# to turn warnings back on, use : options(warn=0)\n",
    "subset_t <- t(subset_without_outlier)\n",
    "prcomp_result <- prcomp(x = subset_t)\n",
    "result_pca <- prcomp_result$x\n",
    "#result_pca\n",
    "rownames_pca_result <- rownames(result_pca)\n",
    "colors <- c()\n",
    "for (row in rownames_pca_result){ colors <- c(colors, ifelse(grepl(\"Case\", row), \"red\", \"blue\")) }\n",
    "plot(result_pca[,1], result_pca[,2], col=colors)\n",
    "legend(\"topright\", c(\"Case\",\"Control\"), pch = 1, col=c(\"red\", \"blue\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise \n",
    "\n",
    "Generate a 3D plot using the function **plot3D::scatter3D()**\n",
    "\n",
    "Tips: \n",
    "- see value of *result_pca* \n",
    "- see ?plot3D::scatter3D\n",
    "- try different phi and theta angles to improve the 3D PCA visualization. Try also 0, 60."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#3D\n",
    "# from packages(\"plot3D\")\n",
    "# Enable this line to see documentation:\n",
    "#?plot3D::scatter3D\n",
    "\n"
   ]
  }
 ],
 "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.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}