2. Data wrangling - Prepping raw clinical and transcript counts data to generate gene level data for DEseq2 using tximport.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Wrangling\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Much of the unique code that you will write for RNAseq analyses and many other data science projects is to extract and reformat data. Let's review our question and pipeline to determine what data we need to extract.\n",
    "\n",
    "## Let's review the question that we will answer:\n",
    "### What are the differentially expressed genes between HER2 positive and Triple-Negative Breast Cancer (TNBC)?\n",
    "Breast cancer is a classic example of the use of precision treatment. Patients with cancerous cells that express the estrogen receptor (ER+) receive hormone therapy and a recent [report](https://www.ncbi.nlm.nih.gov/pubmed/29860917?dopt=Abstract) demonstrates that a subset of those patients based entirely on gene expression do not benefit from chemotherapy. Likewise, patients with cancerous cells with amplified ERBB2 (HER2) are treated with a targeted therapy that has dramatically improved the prognosis of those patients. However, patients with breast cancer cells that do not express ER, HER2, or the progesteron receptor (PR) have what is known as triple-negative breast cancer (TNBC). There are currently no targeted therapies for TNBC breast cancer. One might be interested in what genes are expressed differently in HER2+ vs. TNBC given that HER2+ are typically also ER and PR negative. We will ask that question using the TCGA counts data. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our pipeline\n",
    "In the pipeline image, the red text indicates the names of the R packages we will use to perform the analysis, the yellow boxes indicate the data (raw or processed), and the grey indicates what we will do in each step. Note that we will use [DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html) to identify differentially expressed genes (DEGs). This is a highly cited package that uses a negative binomial distribution to calculate DEGs. However, there are multiple pipelines and programs to do this. [Knoweng](http://education.knoweng.org/sequenceng/) has a nice interface to explore other options, including the strengths and weaknesses of each.\n",
    "![](https://github.com/davidnboone/test-image/blob/master/overview-image.jpg?raw=true \"pipeline overview\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To perform this analysis we need to:\n",
    "1. Extract and calculate HER2 and TNBC status of each sample.\n",
    "2. Extract count and abundance data from those samples.\n",
    "3. The count and abundance data is transformed and currently at the transcript level, so we will need to untransform and collapse to gene level to perform the differential expressional analysis. In this notebook we will prep the data into the proper form to collapse into gene level.\n",
    "\n",
    "___\n",
    "___\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We still need the proper R packages loaded. We will need to start by reloading them with the same code as the previous Notebook.\n",
    "Of course you can use the typical library() command to load each as well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#you can very easily use \"install.packages\" or \"biocLite\" to install the packages and \"library\" to load them\n",
    "#however, instead I am installing a package called pacman that will determine if a package is already installed\n",
    "#if it is not it will install it and after will load it\n",
    "#this method is convenient when sharing code that requires the use of others packages\n",
    "\n",
    "if (!require(\"pacman\")) install.packages(\"pacman\")\n",
    "pacman::p_load(BiocManager, R.utils, data.table, tximport, DESeq2, biomaRt, jsonlite, BiocParallel, ggplot2, gplots, RColorBrewer, devtools, pheatmap)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the working directory to the workspace set up in the previous notebook.\n",
    " You can check to deterimine if you are in the correct directory by examining your working directory. If it is not the proper workspace then change to proper directory.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "getwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#only need to run if you are not in the proper working directory.\n",
    "base_dir <- \"C:/Users/dnb14/Documents\"\n",
    "setwd(base_dir)\n",
    "\n",
    "project_name <- \"TCGA_Her2_TNBC_DEGs\"\n",
    "date <- \"2019_03_01\" #use the same date as notebook 1\n",
    "setwd(project_name)\n",
    "\n",
    "getwd()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "___\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Start the analysis. Extract and calculate HER2 and TNBC status of each patient.\n",
    "To do this you must understand the TCGA Barcode. Let's take a look at what each portion of the barcode means."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TCGA Barcode\n",
    " ![TCGA barcode](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWQAAACOCAMAAAAxbYJCAAABFFBMVEX///8AAADqmZmZmZmWlpbu7u7p6en39/fplJT8/PyTk5Oampr7WEH7YU319fWxsbHCx8fd3d321dX2aVe4wcLV1dXVmZS8vLylpaX9uLH9vrj8nJLY2NjMzMz44OD32tp0dHRpaWnwurqqqqr4dWV6enqMjIz1j4W3t7eHh4d4eHhubm5gYGAjIyMbGxszMzM9PT0pKSkRERG0dnZUVFSXY2M7OztKSkrZjo5NMzPh3+q3d3fJg4PLgHm/vsfMyNYzISFhPz86JibklZX2wcEaEhKbZWV+UlLX1OO+vcvHyb1JMC+5urFbOzvRy+Dr6PLY2c3ImZnPzdPT08bl5tn19ey0tLyDVVXWt7ejfX2umJjTnJbVTmXnAAAZjUlEQVR4nO1dC1ucSLqGNFSlYBa6ZoctZg81C8PCgQZtNUaN8WicnExax4zZmcTZs/v//8f5vgL6Zms33ba5bL95HiVYFFVvfbe6omkbbLDBBhtssMEGG2wwD99/86d2+Oa/PnWRvzwAyd+0wYbkJfD9N3/+918Wx7//viG5Pb7/5u/fdVvgvzcktweQTMzFQTYkLwEg2WGLw9mQvASQ5BbJNyQvgw3Jj4D5JJvjf1+CZGrStoX62jBGMpHeLL59Mfaf9iQ7cd4JJu4A50ySdrm0gSmz7kS7Mrb4wzROjY5hx8CEiCXegV8tnp+JEcndsCz9GUKXx2P/aU8y6Ym44GM3aEo0lsq2BV0YplVGuTl+x5tVrTvApG9FRmo5vNdLikAjSZgUqxZ2SDIpuxqDsjHH0UziONTRHEIdpuUucNuo/DIkE83IIFuqMQqmg8KLTA3yxTuEmm08wiKgVk7NjEKRGRYbCm6KJKCaehPUav77SAFqxgwLfmjMTuHHqmUakiwL4XEoi58kgudJmBaEF0noIskyT4xKNpYgOQmCxGN2aDCeW6HtFlFoULjjhxaN8rB4YMMhI6U2bpLYmhfmkc3SqCgET5M80HIjL+bm4ERg3qiRuJmnUTuJ3ZW1bkgydYFZn7phYFm8dAthCR65MuSpSxJplxUX7UmGlopc1u3JxOP7tsgJibomjQIZykRGPa8n5ufRBqmNP73QS0vNK1OZEJRkYmVZKbQikeHcHFjpwU8zs/JCVL/cFcs05vhAk6MgEZrt85z2SBaT0GSJl7penvdEpTNLkJxIUFwr03yX9wiDtoI60IjbrhbHUO+HJrlQ+cVpkkhNloSERAt6jpmkBTizUGh8XgaaFkkkAy7cEL20lpUrlqkhGdwRo7JEkn2bpywPhEFCbqIke2VMnWXNBU/wScNllqwusQ5Ash+zjog8LclWrMEUUts0c+GWXTDHMqEEnG7QM83SAvuvhd1FskCShU1ML09lB38lK5ZpSHIe2nkUs7jn54kEkqXMSZGmqZO7pmX5SVW8ZUn2ojQBh61ILpHkIIA7BIQ691eswa33hZFBeGLbPY4kg4ENQCPj3IcwqViI5J7EAveSKAe3FMKvBaT/XoyiC+EKcHyOzHhGuGZLkpEwFqYWmBh5elV40Z5kypWhgWw1swpfBYRw3NE8uCMdTT50MEc8yFwjXQHxl8covI3Ae6mXSRNfvQBUgU3PI2z4azXc1ePDAJ6Ht8q06VYvg/u61TzZkPwguI9kJ7jVU9qQvAw2o3CPgMcm2RHxGseG1gTGyUq9/3tIhqzN6Xsrk8wNa9X+UzswwleODqTtu90VsrmHZMg6FlNjkquTbBmPSzLxrZWDRNcyDMOO5bID4/eQHGPWHduXYwm+PJI9y1i5427GPjBhQE5kKZ6B5LueM93YtjudjmHEnlPrCv3SSGbCWJ1kjRLp2gZSkS1jNb7/5s9/vQvf/+Xof2ocfV/fW3VxyycgufMgfUqHSNvCSZOsfSf7sZdpPTbJNDPsYH6yxeDFIM4gzW1Djb/+6c8t8bgk41gKeJ6lh2ic2LBXHd8ZA8+QZlu2o5l8twiyTl5ERVGk8XffrVbktiTLyAp9VjYa33oC1vQN+1YgugJoEKNtdluVg5IF4cV5nuaWy1eb4G9LcjfUZEiihmQrvjf1bZgQHz3skgQq0TT73oNm2oBxbMPOih22tiSLnkZCE0gmvs3NHJSJO76x+CQ9sTp2+1LeD5p1jI7V2jIvCO6j3e+uIhmtSc59MBeRx8tOGmki6XR5YnfKhS2Atw5Hyzwgwoof0gyNwemi3XdXyL0tyVnoC6pFnm1pWa5pvq3JpJsbC7fzg4TJt2H6VsdY2ygMQWFeIfe2JLsW/ow8KwOhBpINTRRWi/5tbKzeq54FR6BhfsC4ZQIQeELuS4eerUlO8Wfp+TnxCyA51WQpAn+huTqE3bHWxIS0oaO9LpY1D0yGtaxzbR3CKWW3CbES4Flz4WGZJwtHZdTqWCsPws0GCyAMsNc2bsvRMC/J8vI9Pkcz6+UQbHG/ziG4WBPJIG3Asr82ltEwL9mPeuRutYRO8PpyD8Bi+GuKMYDleFmWH5lk1zDWt34UZBn909o0BWTZsJeJxh+XZAbFXJtvQnTB5i/shFuDQwc+W6INH5dkYhvrc00IjLUecgBqCgFE40to4uOS7EHPaU293xpmbBjrc37gU7ANGfFaxRmrkNxmW0IF6O+tT5krcHB+6/OtFJyKS3nHMto0ZDuSJ1l937YLBGL2cCP2dwHM8hrtPgHXGoBN6rQRZUXyYrP07OP738fb7+NVd+Ix8v79nPAJTPIaVbmBvc4IQxmMrrG4aZbCUyR7dmeBMaaPVxlgNAdB5cXpiDH2Hv58k907isG6xmM4AN6CgpYgngSyOjb8W1BYeGpYLjSKv5Cv+HB1enE2uLhpVruy96fXg98aTmn35mIwGJzd3McyhQBuPWPrk6/JDGM9w56OYVkgyIrlBcfEhIVz6h3E/OiPXJ0Ojh16fHZz9d6k2scPv2fHmnNaGz96efPrNfjB68HpPeaQgLV8jI2a6Ps8jQrx0KbZszoVFifZcYfPzFWvj1c3A9UOzuAmy64urzL1/2PxAW8Cx8dVK7Hj7qzm4pnHVHfvccJFCGJ8B/TTeuB8id1I5eJzaFTgWhl7kZjn8vS4uTw+uwCcKc7pWfaegjke/ZWezTA80sLetGl0rMdZ3AhkWDha9OBjUabsxhXLLbLGWASJnmfBPmSDUabO8WBwfF1dX1+4V1fuxfEo6fFvtx+XatoC5Mtdn9cfwpTS8zsdH+u1huyJq1ZxtdhZg8N38MTcCaHgYqKfNtYLuT4Df3g9xh0d3H7cUdM3UHPhrV+Uwe2p1WztmFgc1LMNI23VF0GS54/4e8d3/sm5vp40T7O6zbKxZLa97vCC2M271jYH48Vxq5wZiJc9fzL9eMUBB27XLhbM/5rjC1zRWJPcqus7H863Nf727SR+uD/5t9/+81//+se85LOlsw3MuBbkR+jxOVz6Fc3WwwbLf3t6B36cmfyHUYInr1/PTb46GHpYlK01DkOOoXI1Dx3L/O3pk5m4k+TZyZ+si2TqVh0l8UiHxjgCSRYPG8t87iQ7Lq4KfITBodELrZXW/MzC506yFkCvf/nFIUvACZbeUXIXPnuSGV9xV9hngM+e5MVBZwpgc4oOl91grqUFKQ6aFmWEz7QaXA77RTQYnahlenKY//TpVzXJJy9Pjhq+jl5Okjx+2E9F8gkkH6U+qkgO7tZpKNcwNiAjZTRHImrymTv4CZR8WIw5y8ZJL8wty0pTy/brJX40tqwoV9fZoa7v5KP+DJk1Op8VB9sHYdWLlfn+YXl7MCfId7Z2S7c6JCGE5NXRJ8wt+9s7SZW/k+7lEw8qkk/e/vTm1YufXyvWnr3QT4YkYzGLsd3vSPLR819evfnpedUoJ+f6z4pkb2cnnU0zt3Z2dw/rpSOy1OtxO5JG+1F117TKnfL2WRSOv7+7tVP3vkih37882dXHkKt3HeLlHr5PbKvb241s0lm5ifphfMCuHiimWA4OqyQGXHf3qmtk2aouD7CsvITnJiwcknzyS5XkD2TrOTz6siG5LuZoChNIfn1epT5Hll++qh/7MYZbu7POeiJlXW91Pg2UvRq441VxQ7jr7ajL3ak2Ykn1oDo0yTusmbsTBETvcFffO9w56O/6NbE7fUUCXO7aIunvNQ2Z6vrhdAZyW+8bstPX+6Zm63o/t/d1fXIQg+zre7lArfA0b1fvW9LuY8tlkJ3fzbf1yNEIvFLvTZN89ELX3z579pMi9zlWq5FksTUs5ojkt7r+7tmzc0XuS2zM54pkECRs/FtnPtAIXilEqciNm5bXWAS18EMUGyjtXuSH29PVBunYzwQwnTLNg5Lo8wZ5GSGJ3g+4J9EkwSN7Pg/wDU5ZvdVr1vVzlOrpp3N9CxPFe8BsZy8E1Ze7U/IuKhmGEpdQvG1M7u5Bmn19ByUkR7Uhu/1t/Za5eFYJ48kr/cWTJ3/or4YkB7vQkjzY08cl+Qg4Rrvyi/7m6MmzPUj9c03yrgDO9P2pcRwPpBV+Qfv2Nc3Xdypt1Lw+FtfBEsqdHehqsFLfmnjQ2dL3UbvhPgEpO9ibYy4UoJqNRQBxcxWxDrbtxOFDjlKuKXEA1YrwNznUE+YEmI3Ym5JkEBhkjyagdvu6OjUKhDsUu4p7LUDJZ5wc3CYZhBI1//Vb/c3rJ0cnP6v/IsmgLXFdzBHJ8OdnyOrPKPivT05GJPfBOUPNDiYLD5KodBTk0tQcT6CmaUpIoPx0T08xSMOKH+gHEw9mtWERKIU04NOqOxM+vqW5yhvfBoo9ob6xvl1u61PzDbKqLTJZn+oGtdmdEBm2rVeH0bn6lr9dlyfUS9CBSg739ATZ3dEnV/kCyW/08yc1bRgzvGxIHi9mgx+enuu/qMjiqGL3dc25IhmoysE4TDzR10unrgOKcFBzDnLTjzNoxkZLQAsmJTWvVFCjuq6G0CvBnAPRCCjd1/eb1uZIH0UrIlWEEkCZZH/axHdrHQOS1ZMO+o/JIrHGZAHJVlP0RN8P9X5l7Lf0ECu7oxvTJNeuC83Gy+qXIpmNFXOIH56+098d1ST/cZtk9FZbE0t96sZFktGENSRrWeWbk6o0Hijw/qQnD+vXs9rjLUSybLIX2yPBN/Qt3k0Ot5qwAwSBg8WIJh8dIxklmaDbneQKi6LMApJsN8l7eplUYYWmbTeSPDmzMknySfULbO7T/+1uzTCCSPJ5Q/LzGSRrwdakWjfqMJJk0rCAgUuppNUH7xpNhXBhrbVNzfRp/Z4F0lTd13eHU7w9fd/t68PYC7xEmRe7+qGSVzuKorK0HSxgJR1Av6mZoGT9yuEzUWASCDGb9gZ7s+U28UChR7ner8zKHjrpEclMhOrJfyLJb8fNRUMytNZoJtpLMHUkJyS5NhcvJ0iGNwCrXZV7GGDj9pSEgqAEYyQHO+D7lfxqSrT0evkur14klBqqlE5tjBcluaq6re/UdkOxlvV39y0XArBRNFzFEqJSqF2O0YgSKnB8OZpxPW3yPKjSuxo7UMVFx9eHqECZDnB8ibunAg2snFtRUIknqWLTvf978uQnjCoqx/dkjGRf7w+LqdWxbu+Hp2/1VydjTTKL5JQ5dcxuo/c8VPSBDWNjJLuVyCm3iLFaY5jCOjbGByrxELVwjkeS95Hs1iQrSWayJ6DwpSM9lVOCUrh3CH1DCGrRz+GgGy6fYsiWMiDK0IjGMGhq+RKmwVkgKK5T0b7v1NGFhLYBI29VL626FA3JTFZPfqsCX5TNo5/0X8ZJjitJZjIRuNwfU/vkh6fPavvwDkO4MZLjCXMhOxaekGKqsB/tAANPrY1IZlb1XzSEIDqjLZO8ehFXlCi+evo2b0NyZaxAQqFnxNMtvWS1yYSAZi/QzFLfwTCaTcdZUFKUbSfStzz0ul4db4+jW3k+Cz0pxMmgWjTU9zxoxh1P6UAlLDtTgRDY5JeVvfij7lbAf58om7w9LOZ4CAesvjhSqc5VN7wh2a9IxhByYr0LsIphj1urH69iLAbOCHOFCCYAWfCJJ6eHJiBewiKDFhfmwiSbQy0vwKZCBxDpgGBw37CjLdUAwH5lV3toI8YBXYytMjxUcX2h7/V3t7a2dibX60IL7e2HoNYQEkKkv1X2IHmkBKIfFQdKdZ3kcE/f2hnPG3t874C3c/jxCsh79gt0L1798frJj6i5VTHHmqXq8f10fv5Gf/MSRy6gn/gKPeGPtr7dpaYsp+J+bGr9MCzhvRxHLsABHeSOioL3M2lvg5WEau9ijXanYirw3wdhBOIH5BKs/G45d+zeHDaFqSzctuob1cZuFweFbFUQTTX3VL8p225MFQQkNabKBF0QxCFmIbaGybW0Mu3YwKSx4BMkPzl5oW6r+PetunxxBHFyU8xxwRyNXbx6VpkWANINklxjf2qAsB672MFsqnEUbGUzHHof2Xh+fSpgqsYuttB0dqsEc9eB02w4uua4aVLPWhHrcHt7J1EV4XnjP/18eihTWlGZ4CH4eFZR7LpufGunKveLMqznG6WByatMRBpFVc5U+IhxQ1ONwj0/f3f+XPmzk5+fA16qHl9VzIm+khqFe/b23flbZSSOVGqk+0diKCfcv714g7hFWRgqVOMxvl8tejPd8nBnH0/Eh1JhhVx/2hw4Ii+jXBFrZvhgy2VEbFhRRoKgOXliVPsZ48/m/A9QmGPnyo8lZ6Z5Z+nq8eSjo+FwcoVqFI5Nv7IaT349nVoN2nNPCG/monDHnDX8bZLZg8hjoPeU/AvC3JkRRulYPT/nmZHPF/NIDmLbGpsgv49kyh98We9XgjkkO7g+bGxdyT0kO1nHWu9Wuy8W95OsVj10xr6Tcw/JAhdIrHH78xeMe0kGjm1c6T1aA383yaav1iuv/OGfrwZjq4TvI9nEJVJqSf0wOLibZK9eEbkxGDWuRxHVPSSz2KjXrY62EN5F8tO/VGvW2p6z9knBrls/Mamn9PruZTdsMCLiLpLByHLcfROrxf3DjvWdkvyP4TLitocGrg1z9xZfn7XN0hlf0s5wD8xg1JOYfN316UIL6gLDsCRVS1fnH7FkV3alM+knHw90WiqdY8D9NB+ftn2LczYSIDq4ca/cemsXaAW8bmxzhnN2+XGhLD3RpXjW6wKm1gSpjwWacGONJ6Dcg/Hqa0rKTm9uTgejm5OdKsT1xWXbt7DfhqLMBjdX3CFXN+qOMzhTG7uGa7wG2fsF86RqxNq2LSue06WXFkg92Jc4cN1H2Dp6G/S38W09bHDqXl1eZjdnzYfxBr8CBuP7gK4vsg+tX0OafUbA4u/qvfIU7rCLm+zq6iq7uaiUB6T8sk2Yhedvm97cJ4SfMdxl/snsMTkdGQxnkF0igeZlpqhng4vTmwwk+2JIM7s+y35vH27S7hnuJqLXvzaSyuSZYvzDR83hl9nFwGEO7r5ttfzWW/xgQfNhD99tB6x+fXWNGzGrS09txkTreRm8l2BCz6qdQtfHZ+6CRnMSJtB4fDy4yIZLOM3T62tR6wSVGdiMs4usZd7cXphk3D60zsOY7oeZneEubXr864U7sgMe6PfxzaWJBpnRD1c3F0Dz9eDs1G2lzyPwLAOtuBoTJnl6+ntzzTgozE3r9gOS3UVFX3aMB95t0QZQv4uzs4tT0R2ztVT8djziQGOeAKMBxuNqUb90C+zD5eTpHOZEe7H3l+9bq0gbkon/mLs7br/eE5kQ8sNEcfnN6YRcYaKs+/sDFpNdTma2zKFZLUjG3fyfJLaowZwPH5ypOlJxNSlYzPxAphOthunYsD1akSwfZ6d5K5D2gdqjA0he/NQ4PLX0yxm7+HzQIrrQ1KEcj/sZna8D0JdrEZbhCSWbOai2gNi3zYFsGCqv+STCrxBtv+6B56t+wlj5y4Tp2+2O65c4B7Wmwny1MHk77WeZ9XXO9JEH/2T0CqDx+s4JXKwAXX/iRBk2vVx0uWzzoje5TRF3/Mb1ZgRrxjQ9X2cwS7MFZlLWCNIz0nysgjyc6EyxbLlJhUQ4HWPyBi5Nq65nLU+z13JaVAPH+6TTfCThTjrmSOiURBXLSUCvq2UpSFAcaMLLXNIpEt+UnubEtqcJIkSGS0g9v0u7QtiEGr3w0w38rh0kISyXmbATQt3U02jOIX7PM6qJXkpkERXpMrIc+nEhnTyxCqeM0ki6hZ3RPNb8xM9pIf0yDH0ti+yiG5dJYjA3Sb+gmfu2IGEnsRyr9HsiTuyCmqFnJrFdkqAQbteUhTt/omcGohyMoCi0IKGlyyCuTaFvm2c8RCYjGVuaVN8Pi7w4Mru4OWWZsyiX+3TrEM5jza+Swu8CBbYmSSK1MKA96eVZYtP663ThckeKq4/cuZYWgNx66D3TGO2yVGvxS0+RTIDx0PNT/GSbFqdzckTgGk3mEtKYsLztVP+kwHSHGwjbH7beClxtok/R3/e4FnJkIsLOfq/e3bjcd7pK9HGiZ4qCRsrdWeDYksyLPB7TUropt3MnkqQQ4PFA4rX4/lMNKiQxfg6Sd5u0ZUuHwYuJVsmGe0nkGo/bxvf2kOQcB6osy0+QCdKLpSX8XNpQp3C5VZEpejFiFXkuq8NJbMjHBwsSFjazuBuFYPuzMDR8P9a8EKR+kc8HxlC+3FenzlBcth+1G48309LFtf+met7B78PhWLnDzCwJqLq1HphqB0aMIkFsi5saeCduW7FJ4xT9/WqfkjSZqVXL/53afFbbCGJLzQQ4jvorRjTOIhUMCkJAGCHe9sCnkrYkizxF25SXMSN2kntZTlOpmTlJo6IQPE3StQY41XQFYzQlqNfVBhFn9JeHhr+IBZ6BSIoQF2CZPdcreEuSnVBkoSZKXyQ0TXkvEznoBuuGxBQJJ6nIynV8w/AWWCePHqNrJNp+87VGbhvwZMf3EscMW5JMjTIsIlOUmtcz4cmkCyTzUKChx/ySNHykIVEn/qynEWQUgWkxgGQTelHtSOaRNFmvC4bYC5HkXODnJLMI/D3Q7pilZdKvYmfTqgAmKHbBSWJ3CtGOZDAV+B1UIDmInDT3ExsjGghh4TbQzvzcz+ctsPuPABMYgxGicSFloHlt1JtgGGVywjUzo6QruMBDcHih4qAuBBee+FzWMX9t8BeJzzdYCU667o9kbaAxvrHCG2ywwQYbbLDBBhtssMEGG2ywwX8M/h/nvaBbiW8b9gAAAABJRU5ErkJggg== \"barcode\")\n",
    "For our purposes we will focus on the first 12 alphanumeric digits of this code. \"TCGA-XX-XXXX-XX\". The first 6 Xs provide an anonymous identifier unique to each patient. The final 2 Xs (highlighed in a red box) define what type of sample it is, which is defined [here](https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes). For the breast cancer data that we downloaded there are only three codes, \"-01\" is a solid tumor, \"-06\" is a metastatic tumor, and \"-11\" is normal solid tissue. For asking what are the difference between TNBC and HER2+ primary cancers, we will only be interested in the samples with \"-01\" as these are primary tumors. This is our first piece of information to extract data. However, we also need to know if the tumors are TNBC or HER2+. We can infer this based on the clinical data that we downloaded."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining HER2+ and TNBC tumors\n",
    "Read in the clinical data and examine it by running the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################################################################\n",
    "###################  Define HER2 positive and TNBC BC tumors ###############################\n",
    "############################################################################################\n",
    "\n",
    "\n",
    "# Read clinical data. Note from experience I know that there are both NAs and blank spaces. The na.strings = c(\"\", \"NA\") command will force both to be NAs in the dataframe\n",
    "brca.clin.data <- read.delim(\"./raw_data/tcga_brca.clinical\", header = TRUE, stringsAsFactors = FALSE, na.strings = c(\"\",\"NA\"))\n",
    "head(brca.clin.data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table is set up so that patients (notice the TCGA barcode) are listed as rows and clinical features as columns. To see a list of all clinical features by examing the column names. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(brca.clin.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Where are HER2 and TNBC statuses?\n",
    "To determine TNBC we need to know if each tumor is positive or negative for estrogen receptor (ER), progesterone receptor (PR), and HER2. ER is determined by immunohistochemistry (IHC) and is directly reported for all patients under \"breast_carcinoma_estrogen_receptor_status.\" PR is determined by immunohistochemistry and is directly reported for all patients under \"breast_carcinoma_progesterone_receptor_status.\" However, HER2 is not directly reported for each patient. Instead, HER2 is determined clinically by two different tests (IHC and in situ hybridization) that are reported (if available) under \"lab_proc_her2_neu_immunohistochemistry_receptor_status\" and \"lab_procedure_her2_neu_in_situ_hybrid_outcome_type.\" Often IHC is performed first and if an indeterminable result is obtained then in situ hybridization is performed. However, if you examine the data closely you will notice that in some instances IHC and in situ hybridization will be discordant. We will throw out those samples for our analysis as we cannot properly determine HER2 status. Otherwise if one of the two tests has a \"positive\" or \"negative\" result and the other test is NA or indeterminable, we will use the single test result. \n",
    "\n",
    "The clinical data above has more data than we are interested in, which can make it a bit more difficult to work with. Run the code below to select only the columns that we need to determine HER2 and TNBC status. These are the columns that provide information on ER, PR, and HER2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#code to pull out ER and PR final statuses and HER2 tests used to calculuate HER2 final status to make the clinical data easier to work with\n",
    "brca.clin.selected <- brca.clin.data[, c(\"sampleID\", \"breast_carcinoma_estrogen_receptor_status\", \"breast_carcinoma_progesterone_receptor_status\", \"lab_proc_her2_neu_immunohistochemistry_receptor_status\", \"lab_procedure_her2_neu_in_situ_hybrid_outcome_type\")]\n",
    "head(brca.clin.selected)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### pseudocode to determine HER2 status\n",
    "In the psuedocode below I abbreviate the column names lab_proc_her2_neu_immunohistochemistry_receptor_status as HER2 IHC and lab_procedure_her2_neu_in_situ_hybrid_outcome_type as HER2 FISH. \n",
    "\n",
    "Note that we replace NAs to \"unknown\" to make execution easier.\n",
    "\n",
    "We will define each sample as HER2 positive, negative, or unknown. Unknown means that either the FISH and IHC are discordant or there was not a positive or negative call for either.\n",
    " \n",
    "\n",
    "\n",
    "- if HER2 IHC is \"positive\" AND HER2 FISH is NOT \"negative\" then the sample is HER2 positive, else\n",
    "- if HER2 FISH is \"positive\" AND HER2 IHC is NOT \"negative\" then the sample is HER2 positive, else\n",
    "- if HER2 IHC is \"negative\" AND HER2 FISH is NOT \"positive\" then the sample is HER2 negative, else\n",
    "- if HER2 FISH is \"negative\" AND HER2 IHC is NOT \"positive\" then the sample is HER2 negative, else\n",
    "- we don't know\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set NAs to \"unknown\" to allow for easier logic in following steps\n",
    "brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status[is.na(brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status)] <- \"unknown\"\n",
    "brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type[is.na(brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type)] <- \"unknown\"\n",
    "\n",
    "#create a new column that defines HER2 status as \"Positive\", \"Negative\", or \"unknown\"\n",
    "brca.clin.selected$HER2 <- ifelse(brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status == \"Positive\" &\n",
    "                                    brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type != \"Negative\", \"Positive\",\n",
    "                                  ifelse(brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type == \"Positive\" &\n",
    "                                           brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status != \"Negative\", \"Positive\",\n",
    "                                         ifelse(brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status == \"Negative\" &\n",
    "                                                  brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type != \"Positive\", \"Negative\",\n",
    "                                                ifelse(brca.clin.selected$lab_procedure_her2_neu_in_situ_hybrid_outcome_type == \"Negative\" &\n",
    "                                                         brca.clin.selected$lab_proc_her2_neu_immunohistochemistry_receptor_status != \"Positive\", \"Negative\", \"unknown\"))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the logic\n",
    "Let's check that it worked in two different ways: \n",
    "1. spot checking the logic by looking directly at the table.\n",
    "2. calculating the percent of samples that are HER2+, which should be ~15-20%. This is the known percentage of women diagnosed with HER2+ breast cancer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(brca.clin.selected)\n",
    "\n",
    "#calculate number of HER2 pos, neg, and unknowns\n",
    "num.her2.neg <- nrow(brca.clin.selected[brca.clin.selected$HER2 == \"Negative\", ])\n",
    "num.her2.pos <- nrow(brca.clin.selected[brca.clin.selected$HER2 == \"Positive\", ])\n",
    "num.her2.unk <- nrow(brca.clin.selected[brca.clin.selected$HER2 == \"unknown\", ])\n",
    "\n",
    "#calculate percentage of HER2+ tumors\n",
    "percent.her2.pos <- 100* num.her2.pos / (num.her2.neg + num.her2.pos)\n",
    "\n",
    "message(\"number of HER2 negative samples\")\n",
    "num.her2.neg\n",
    "\n",
    "message(\"number of HER2 positive samples\")\n",
    "num.her2.pos\n",
    "\n",
    "message(\"number with indeterminate HER2 status\")\n",
    "num.her2.unk\n",
    "\n",
    "message(\"percent of HER2 positive samples\")\n",
    "percent.her2.pos\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We cannot use all of the samples.\n",
    "From the code above, we can see that there are 857 definitively HER2 negative tumors, 206 HER2 positive tumors, and 184 that we do not know. For the purposes of this exercise, we will exclude any unknowns. This makes our dataset smaller, but cleaner. We may even lose more samples after determining if samples are triple negative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now you write the logic to determine if each sample is TNBC or not.\n",
    "Note, if a sample is TNBC then ER, PR, **AND**, HER2 must all be negative. If any of them are positive, then they are not TNBC. There also will be some that we cannot determine because the known ER, PR, and HER2 will be negative, but at least one of the three is unknown. Store the values in a column named \"TNBC\" in the brca.clin.selected data frame. Have \"NOT_TNBC\" indicicate that the sample is NOT_TNBC, \"TNBC\" indicate that the sample is TNBC, and \"unknown\" mean it is unclear with the given data.\n",
    "\n",
    "**Hint** be sure to change ER and PR NAs to \"unknown\" to make the logic easier like in the code above. The results should indicate that 176 samples are TNBC and 973 are not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brca.clin.selected$breast_carcinoma_estrogen_receptor_status[is.na(brca.clin.selected$breast_carcinoma_estrogen_receptor_status)] <- \"unknown\"\n",
    "brca.clin.selected$breast_carcinoma_progesterone_receptor_status[is.na(brca.clin.selected$breast_carcinoma_progesterone_receptor_status)] <- \"unknown\"\n",
    "\n",
    "#if sample is TNBC mark it as \"TNBC\", if it is not mark it as \"NOT_TNBC\", if it is unknown mark it as \"Unknown\"\n",
    "brca.clin.selected$TNBC <- \n",
    "\n",
    "\n",
    "#determine number of TNBC and NOT_TNBC\n",
    "nrow(brca.clin.selected[brca.clin.selected$TNBC == \"TNBC\", ])\n",
    "nrow(brca.clin.selected[brca.clin.selected$TNBC == \"NOT_TNBC\", ])\n",
    "\n",
    "#we are interested in the number of TNBC and number of HER2+\n",
    "#count number of HER2+ tumors\n",
    "nrow(brca.clin.selected[brca.clin.selected$HER2 == \"Positive\", ])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the clinical data table\n",
    "Now that we extracted and computed all of the appropriate data from the larger data table, we should save the new table with the generated feature into the dated directory. That way we can do not need to run the above code again, but can simply read in the table if we need it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write clinical table with HER2 and TNBC features in analysis directory (date) so we can just load later instead of running the above code\n",
    "write.table(brca.clin.selected, file = file.path(date, \"brca.clin.data.her2.tnbc.txt\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if you need to read the table you just saved, use this code, just remove the comment.\n",
    "#do not do this unless you need to read the table\n",
    "#brca.clin.selected <- read.table(file = file.path(date,\"brca.clin.data.her2.tnbc.txt\"), header = TRUE, stringsAsFactors = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "___\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Read TCGA count and abundance data and parse for only HER2+ and TNBC Breast Cancer data\n",
    "The downloaded data contains counts from ~197,000 transcripts across ~10,500 samples resulting in a nearly 15GB file. Simply reading this file using base R would require more RAM than your personal computer likely has, so instead we will use the fread function with the select command in the data.table package to only read the breast cancer data (both counts and tpm) that we are interested in exploring. fread is also much faster than read.table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################################################################################################\n",
    "######################### read counts data and parse only BRCA data ###############\n",
    "###################################################################################################################\n",
    "\n",
    "\n",
    "#note fread with selected columns is required because total tcga count file is too large to read into RAM\n",
    "#the columns selected are those of breast cancer samples with both clinical and counts data\n",
    "#fread can also not handle duplicates so those must be removed see the dups<- line\n",
    "\n",
    "#this line reads the first 10 rows of the entire tcga counts data. This is used to find BRCA sample names and duplicates for fread\n",
    "head.tcga.counts <- fread(\"./raw_data/tcga_countdata\", nrows = 10)\n",
    "\n",
    "#lets look at the first 10 columns of these 10 rows\n",
    "head.tcga.counts[1:10, 1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the count data, note that each column is a patient sample and each row is a transcript with an [Ensembl ID](https://useast.ensembl.org/Homo_sapiens/Info/Index). Note Ensembl transcripts begin with ENST and Ensembl gene IDs begin with ENSG. We want all of the transcripts, but only need IDs corresponding to the 176 TNBC samples and 206 HER2+ samples. To extract these we need to use the gene IDs corresponding to the TNBC and HER2+ samples from the clinical datatable. Let's look at how the format of those IDs correspond to those in the counts table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lets look at these 10 rows\n",
    "head.tcga.counts[1:10,1:10]\n",
    "\n",
    "#and compare the BRCA TCGA IDs from the clinical data\n",
    "head(brca.clin.selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We need to intersect the columns of the entire counts file with those of the BRCA clinical data to provide a list of IDs that is only from HER2+ and TNBC breast cancer\n",
    "**note** If we compare the IDs we see they are in the same format with all uppercase letters and a hyphen seperating the sections of the ID. This is not always the case. Sometimes the IDs are writen tcga.xx.xxxx.xx. If that was the case, we would need to make them the same format before running the intersection command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the intersection of the clinical data samples and counts data is used because we only want counts data\n",
    "#from samples with clinical data.\n",
    "#note we add \"sample\" to the front because it is not in the names of brca.clin.selected but it is the column name \n",
    "#for the transcripts in tcga.counts (see above) \n",
    "#and we need that column to get transcript name in our final count file\n",
    "\n",
    "#note we will only read data from samples that are either HER2+ or TNBC\n",
    "#this command will find all samples that are HER2+ or TNBC AND have counts data\n",
    "#it intersects theh sample names from the clinical datatable with the column names of the counts datatable\n",
    "brca.counts.tnbc.her2.only <- c(\"sample\", intersect(brca.clin.selected[brca.clin.selected$HER2 == \"Positive\" | brca.clin.selected$TNBC == \"TNBC\", 1], colnames(head.tcga.counts)))\n",
    "\n",
    "#this list contains tumors and normals. For this analysis, we will focus on tumors only\n",
    "#Remember, tcga IDs define tumors with \"-01\" at the end of the ID and \"-11\" for normals\n",
    "#we can use grep and the $ regexp to pull out only the tumors\n",
    "#but don't forget that we need \"sample\" too!\n",
    "brca.counts.tnbc.her2.tumors <- c(\"sample\", brca.counts.tnbc.her2.only[grep(\"-01$\", brca.counts.tnbc.her2.only)])\n",
    "\n",
    "#examine the names and calculate how many we have now\n",
    "brca.counts.tnbc.her2.tumors\n",
    "length(brca.counts.tnbc.her2.tumors)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lost more data\n",
    "Note that 339 is less than the sum of the HER2+ (206) and TNBC (176) patients calculated from the clinical data. The loss is because some of the clinical data referred to patients without a primary tumor with RNAseq data. Remember TCGA data contains multiple types of omic data (mutation, rnaseq, proteomics, etc), so this loss is not uncommon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing duplicates before fread\n",
    "In the TCGA data a few samples had more than one section sequenced. This is indicated in the portion part of the TCGA ID code (see above). However, when reanlyzed the TCGA IDs were truncated, so the IDs are now duplicated when multiple sections of the same sample were sequenced. This will cause an error when using fread with the select command to choose the brca samples. Accordingly, we will proceed with the duplicates removed. Other functions can work with duplicates, but for the purpose of this exercise we will stick with the much faster fread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fread can also not handle duplicates so those must be removed\n",
    "\n",
    "#calculates the duplicate samples\n",
    "dups <- colnames(head.tcga.counts)[duplicated(colnames(head.tcga.counts))]\n",
    "\n",
    "#removes the duplicates from the list of IDs of TNBC OR HER2+ tumors with rnaseq count data\n",
    "brca.samples.final <- setdiff(brca.counts.tnbc.her2.tumors, dups)\n",
    "\n",
    "#calculate the number of duplicates\n",
    "number_duplicates <- length(brca.counts.tnbc.her2.tumors) - length(brca.samples.final)\n",
    "number_duplicates\n",
    "\n",
    "#print to screen the final list of HER2+ OR TNBC tumors that are not duplicates and have rnaseq\n",
    "#count data\n",
    "brca.samples.final\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** that for this dataset there were no duplicates, but if you use this data to answer another questions, you might encounter this issue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally lets read the counts and tpm breast cancer and save it. We will save immediately, so that we will not need to run these memory and time intensive lines again. This may take a few minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this is the command that reads in all counts data of BRCA non-duplicated samples with clinical data\n",
    "#it also converts the data.table into a dataframe to use with base R functions\n",
    "\n",
    "\n",
    "brca.counts.df <- data.frame(fread(\"./raw_data/tcga_countdata\", select = brca.samples.final), row.names = 1)\n",
    "save(brca.counts.df, file = file.path(date, \"brca.her2.known.counts.Rda\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We need to do the same for the abundance (tpm) data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "brca.tpm.df <- data.frame(fread(\"./raw_data/tcga_tpmdata\", select = brca.samples.final), row.names = 1)\n",
    "save(brca.tpm.df, file = file.path(date, \"brca.her2.known.tpm.Rda\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "___\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. [DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html) requires raw counts at the gene level. \n",
    "This means that we need to:\n",
    "- a) determine if the counts that we downloaded have been transformed, and if they have untransform them.\n",
    "- b) collapse the raw and untransformed transcript count data to gene level by using tximport"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3a. Undo the transfomation, if any, of the count and abundance data\n",
    "To determine if the count data was transformed, examine the metadata. The metadata is java script, so we will use the fromJSON() function in R to read these files.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################################\n",
    "############### DESeq2 only works on non transformed values need to look at metadata ##########\n",
    "###############################################################################################\n",
    "\n",
    "BRCA.count.metadata <- fromJSON(\"./metadata/tcga_count_metadata\")\n",
    "BRCA.count.metadata\n",
    "\n",
    "BRCA.tpm.metadata <- fromJSON(\"./metadata/tcga_tpm_metadata\")\n",
    "BRCA.tpm.metadata\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examining the $unit value from reading the metadata above, we can see that the counts are actually log2(count + 1) and the abundance (tpm) data are log2(tpm+0.001). It is common when log tranforming data to add a small number to each value before calculating the log to prevent an error when calculating the log of 0. log2 tranformation is the most common for expression data as raw values from gene to gene vary by quite a bit naturally, but not enough for a log10 transformation.\n",
    "\n",
    "We need to undo this transformation for the counts and abundance data, check the new table to be sure the values were untransformed, and then save the resulting dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#under unit it describes the transformation as log2(counts+1)\n",
    "#need to undo this transformation\n",
    "brca.untran.counts.df <- 2 ^ brca.counts.df - 1\n",
    "brca.untran.tpm.df <- 2 ^ brca.tpm.df -.001\n",
    "\n",
    "\n",
    "\n",
    "#lets look at values before and after transformation to make sure\n",
    "message(\"counts raw\")\n",
    "brca.counts.df[1:5,1:5]\n",
    "\n",
    "message(\"counts untransformed\")\n",
    "brca.untran.counts.df[1:5,1:5]\n",
    "\n",
    "message(\"tpm raw\")\n",
    "brca.tpm.df[1:5,1:5]\n",
    "\n",
    "message(\"tpm untransformed\")\n",
    "brca.untran.tpm.df[1:5,1:5]\n",
    "\n",
    "#save the tables\n",
    "save(brca.untran.counts.df, file = file.path(date, \"df.unlogcounts_her2_tnbc.Rda\"))\n",
    "save(brca.untran.tpm.df, file = file.path(date, \"df.unlogtpm_her2_tnbc.Rda\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3b. Prep the transcript-level data to collapse to gene level using [Tximport](https://bioconductor.org/packages/release/bioc/html/tximport.html)\n",
    "DESeq2 calculates differential *gene* expression not *transcript* expression. There are packages that will identify differentially expressed exons or transcripts. However, these are still largely unreliable because many isoforms share the same exons, making accurate mapping difficult. Accordingly, for DESeq2, we need to collapse the transcript counts to gene counts. Tximport is a package that can do this for us. I highly recommend reading the [manual](https://bioconductor.org/packages/release/bioc/html/tximport.html) before starting. \n",
    "\n",
    "Tximport requires:\n",
    "1. counts or tpm data (or output from pseudoaligners Salmon or Kallisto)\n",
    "2. knowledge file that maps each transcript to gene (found in a transcriptome annotation file). \n",
    "\n",
    "The counts and tpm data are ready. However, we will need to create the knowledge file with the appropriate versions of the genome, transcriptome, and gene name. There are many different gene naming conventions. Also, annotation of the human genome and transcriptome is fluid and changes as more discoveries are made. Hence, it is critical when you perform analyses to detail the version of the genome and annotation file of the transcriptome. The authors did that for all users and is published in the metadata that we already downloaded. To determine what the authors used for reanalyzing this data, let's explore their metadata again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################################################################\n",
    "########### Convert transcript to gene level for counts data  ##################################\n",
    "################################################################################################\n",
    "\n",
    "#because transcript level data is still unreliable\n",
    "#find annotation files used from the public data\n",
    "#this is in the metadata for counts and tpm\n",
    "BRCA.count.metadata\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### create tximport knowledge file from .gtf used by the authors\n",
    "There are other ways to do this, but we will focus on a method that will familiarize you with a .gtf - the main type of genome or transcriptome annotation file. From the $probeMap value above, you can see that the authors used gencode v.23 annotation. [Gencode publishes the release history of the annotation](https://www.gencodegenes.org/human/releases.html). If you go to the link, you will notice that there are 2 Ensembl versions of gencodev.23. You can run tximport with both of them and the wrong version will display an error of unused transcripts. However, from experience, I know that these authors used the Ensemblv.81 version. You can also find hints from the publication itself. Let's download Ensebl v.81 and examine it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#download the file from the server\n",
    "download.file(url = \"ftp://ftp.ensembl.org/pub/release-81/gtf/homo_sapiens/Homo_sapiens.GRCh38.81.gtf.gz\", destfile = \"./raw_data/ensembl81.gtf.gz\", method = \"libcurl\")\n",
    "\n",
    "#decompress\n",
    "gunzip(\"./raw_data/ensembl81.gtf.gz\")\n",
    "\n",
    "#read as a data table\n",
    "df81.tmp <- fread(\"./raw_data/ensembl81.gtf\", data.table=options(datatable.fread.datatable=FALSE))\n",
    "head(df81.tmp, nrows = 10)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class(df81.tmp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find information about .gtf files [here](https://www.gencodegenes.org/pages/faq.html). These annotation files contain all of the information about the transcriptome structure. The first column is the chromosome location. The second is the annotation source of the gene, transcript, or exon. The 3rd column indicates if the structure being described is a gene, transcript, or exon (more on this below). The 4th and 5th are the start and stop stop respectively. The 7th column indicates which strand the transcript is on. Finally, the 9th column is a semicolon seperated group of nomenclature including the Ensembl gene and transcript ID. \n",
    "\n",
    "Back to the third column, note the structure from row to row. Information about a gene is first given (see line 1 from the output above for gene ENSG00000223972) including the start and stop position of the entire gene, but ignores transcript and exon/intron structure. The next row (see line 2 from above; transcript ENST00000456328 of gene ENSG00000223972) will then be the first transcript of the preceeding gene again with the start and stop positions as well as the IDs for each, again while ignoring the exon/intron structure. However, the next set of rows will describe the exon/intron structure. In the above example rows 3, 4, and 5 indicate the 3 exons making up the first transcript (ENST00000456328). After that structure is complete the next line indicates the second transcript (ENST00000450305) of the gene (ENSG00000223972). We are interested in creating a table with two columns- transcript and gene. So, we will start by extracting all transcripts. We can do this by by extracting all rows with the third column == to transcript.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the following will rename the columns appropriately\n",
    "colnames(df81.tmp)<- c(\"chr\",\"source\", \"level\", \"start\", \"end\", \"xxx\", \"strand\", \"xx\", \"other\")\n",
    "\n",
    "#in the third column it specifies if the element being described is a gene, transcript, or exon. We want to map all transcripts to gene names\n",
    "#so we need to extract transcripts with the line of code below\n",
    "df81.transcripts <- df81.tmp[df81.tmp$level==\"transcript\", ]\n",
    "\n",
    "head(df81.transcripts)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As you see, we now have a dataframe with only transcripts. Next, notice that the last column has a lot of information including gene and transcript names delimited by a \";\". We can split that column and use regular expressions to pull out all gene names, which start with \"ENSG\" and all transcript names, which start with \"ENST\". We will look at the dataframe after each call to see the changes. **note** when running this code, you will receive the below warning message that tells us that the pieces of information in the last column are not the same length. However, you can ignore this because we will only be extracting ENSG and ENST IDs, which are in every row. The error:\n",
    "```\n",
    "Warning message in .Method(..., deparse.level = deparse.level):\n",
    "“number of columns of result is not a multiple of vector length (arg 1)”\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the gene and transcript names are in the final column, which is itself delimited by a semicolon.\n",
    "#split the column based on ';'\n",
    "transandgene81.df <- data.frame(do.call('rbind', strsplit(as.character(df81.transcripts$other),';',fixed=TRUE)))\n",
    "\n",
    "head(transandgene81.df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice all of the information that was semicolon delimited in the \"other\" column are now columns of this new data frame of all of the transcripts in this version of the transcriptome. It is important to remember that this annotation is evolving, which is why there are different versions of the transcripts and genes. This information is sometimes contained within the IDs themselves in the form ENSGXXXXXXXXXXX.Y where Xs are numbers specifying the specific gene and the Y, which follows a \".\" specifying the annotation version of that gene. For now we will ignore versions.\n",
    "\n",
    "Now run the code below that uses regular expressions to pull out the gene and transcript IDs. The regular expression we are using is \"ENSG.\\*\" which finds anything in the column that starts ENSG and then is followed by any number of characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use regexpressions to pull out the names of the transcript and genes. All genes start with \"ENSG\" all transcripts with \"ENST\"\n",
    "transandgene81.df$gene <- regmatches(transandgene81.df[,1], regexpr(\"ENSG.*\", transandgene81.df[,1]))\n",
    "transandgene81.df$transcript <- regmatches(transandgene81.df[,3], regexpr(\"ENST.*\", transandgene81.df[,3]))\n",
    "\n",
    "#for the tximport knowledge file we only need the transcript and gene names in that order\n",
    "tx2gene81 <- transandgene81.df[, c(\"transcript\", \"gene\")]\n",
    "head(tx2gene81, nrows = 10)\n",
    "\n",
    "#save the tximport knowledge file\n",
    "save(tx2gene81, file = file.path(date, \"tx2gene81.Rda\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We  completed all of the data wrangling necessary to run tximport and DESeq2 in Notebook 3. Similar steps will be followed if anylyzing your own data.\n",
    "\n",
    "1. How wouuld you do the data wrangling, if you were interested in performing a DEG analysis on tumor vs normal adjacent tissue in breast cancer?\n",
    "\n",
    "2. How would you do the data wrangling, if you were interested in another tumor type like sarcomas?\n",
    "\n",
    "___\n",
    "___\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## If you are continuing to Notebook 3 without restarting, you may need to clean up some memory. \n",
    "This is a handy trick to remember anytime you work with large data sets to free up RAM to improve speed and performance. If you are starting fresh after a restart, you will not need to run these steps.\n",
    "\n",
    "If you compiled all of the steps into a single notebook or R script, you can save some time by keeping data in the global environment that you will need for tximport.\n",
    "\n",
    "For tximport we only need \n",
    "1. the unlogged transcript counts data\n",
    "2. the unlogged transcript tpm data\n",
    "3. tx2gene81 \n",
    "\n",
    "After tximport, we will need the clinical data for DESeq2\n",
    "4. brca.clin.selected\n",
    "\n",
    "we can find everything loaded in the global environment with the ls() command, remove all the files except the 3 below with rm() and then clean the ram with gc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if you run this \n",
    "\n",
    "# to view what is in the global environment before clearing\n",
    "message(\"data in global env before clean\")\n",
    "ls()\n",
    "\n",
    "#to remove from the environment everything except the 4 that we need and information about the start date of the analysis to save properly\n",
    "rm(list = setdiff(as.vector(ls()), c(\"brca.clin.selected\", \"brca.untran.counts.df\", \"tx2gene81\", \"brca.untran.tpm.df\", \"date\")))\n",
    "\n",
    "#to make sure it worked\n",
    "message(\"data in global env after clean\")\n",
    "ls()\n",
    "\n",
    "#to free up ram\n",
    "gc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to free up the most ram you can completely clear the global environment by running the commands below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "message(\"data in global env before clean\")\n",
    "\n",
    "#to view what is in the global environment before clearing\n",
    "ls()\n",
    "\n",
    "#to remove everything from the global environment\n",
    "rm(list=ls())\n",
    "gc()\n",
    "\n",
    "#to make sure it worked\n",
    "message(\"data in global env after clean\")\n",
    "ls()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}