2.TCGA-process.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert downloaded TCGA datasets into sample × gene matrices\n",
    "\n",
    "This notebook is updated to include the data from the [TCGA PanCanAtlas April 2018 updates](http://www.cell.com/pb-assets/consortium/pancanceratlas/pancan/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "import os\n",
    "\n",
    "import pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read gene information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th>symbol</th>\n",
       "      <th>description</th>\n",
       "      <th>gene_type</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>A1BG</td>\n",
       "      <td>alpha-1-B glycoprotein</td>\n",
       "      <td>protein-coding</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>A2M</td>\n",
       "      <td>alpha-2-macroglobulin</td>\n",
       "      <td>protein-coding</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               entrez_gene_id symbol             description       gene_type\n",
       "entrez_gene_id                                                              \n",
       "1                           1   A1BG  alpha-1-B glycoprotein  protein-coding\n",
       "2                           2    A2M   alpha-2-macroglobulin  protein-coding"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load genes\n",
    "path = os.path.join('download', 'genes', 'genes.tsv')\n",
    "gene_df = (pandas.read_table(path, dtype='str')\n",
    "    .set_index('entrez_gene_id', drop=False)\n",
    "    [['entrez_gene_id', 'symbol', 'description', 'gene_type']]\n",
    ")\n",
    "gene_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load gene updater\n",
    "path = os.path.join('download', 'genes', 'updater.tsv')\n",
    "updater_df = pandas.read_table(path)\n",
    "old_to_new_entrez = dict(zip(updater_df.old_entrez_gene_id,\n",
    "                             updater_df.new_entrez_gene_id))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>symbol</th>\n",
       "      <th>chromosome</th>\n",
       "      <th>entrez_gene_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(FM-3)</td>\n",
       "      <td>chr2</td>\n",
       "      <td>10316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>(IV)-44</td>\n",
       "      <td>chr14</td>\n",
       "      <td>28337</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    symbol chromosome  entrez_gene_id\n",
       "0   (FM-3)       chr2           10316\n",
       "1  (IV)-44      chr14           28337"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load chromosome-symbol to entrez_gene_id mapping\n",
    "path = os.path.join('download', 'genes', 'chromosome-symbol-mapper.tsv')\n",
    "chr_sym_map_df = pandas.read_table(path)\n",
    "chr_sym_map_df.chromosome = 'chr' + chr_sym_map_df.chromosome\n",
    "chr_sym_map_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read sample information\n",
    "\n",
    "This file contains sample information. See the [online documentation](https://xenabrowser.net/datapages/?dataset=EB%2B%2BAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.xena&host=https%3A%2F%2Fpancanatlas.xenahubs.net) for more information. For more details on curation refer to [Liu et al. 2018](https://doi.org/10.1016/j.cell.2018.02.052 \"An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>acronym</th>\n",
       "      <th>disease</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACC</td>\n",
       "      <td>adrenocortical cancer</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BLCA</td>\n",
       "      <td>bladder urothelial carcinoma</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  acronym                       disease\n",
       "0     ACC         adrenocortical cancer\n",
       "1    BLCA  bladder urothelial carcinoma"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = os.path.join('mapping', 'diseases.tsv')\n",
    "disease_df = pandas.read_table(path)\n",
    "disease_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Code</th>\n",
       "      <th>Definition</th>\n",
       "      <th>Short Letter Code</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>01</td>\n",
       "      <td>Primary Solid Tumor</td>\n",
       "      <td>TP</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>02</td>\n",
       "      <td>Recurrent Solid Tumor</td>\n",
       "      <td>TR</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Code             Definition Short Letter Code\n",
       "0   01    Primary Solid Tumor                TP\n",
       "1   02  Recurrent Solid Tumor                TR"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Data from https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes\n",
    "path = os.path.join('mapping', 'tcga_sampletype_codes.csv')\n",
    "sampletype_codes_df = pandas.read_csv(path, dtype='str')\n",
    "sampletype_codes_dict = dict(zip(sampletype_codes_df.Code, sampletype_codes_df.Definition))\n",
    "sampletype_codes_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = os.path.join('download', 'Survival_SupplementalTable_S1_20171025_xena_sp.tsv.gz')\n",
    "\n",
    "# Mapping to rename and filter columns\n",
    "renamer = collections.OrderedDict([\n",
    "    ('sample', 'sample_id'),\n",
    "    ('_PATIENT', 'patient_id'),\n",
    "    ('cancer type abbreviation', 'acronym'),\n",
    "    ('__placeholder__', 'disease'),\n",
    "    ('age_at_initial_pathologic_diagnosis', 'age_diagnosed'),\n",
    "    ('gender', 'gender'),\n",
    "    ('race', 'race'),\n",
    "    ('ajcc_pathologic_tumor_stage', 'ajcc_stage'),\n",
    "    ('clinical_stage', 'clinical_stage'),\n",
    "    ('histological_type', 'histological_type'),\n",
    "    ('histological_grade', 'histological_grade'),\n",
    "    ('initial_pathologic_dx_year', 'initial_pathologic_dx_year'),\n",
    "    ('menopause_status', 'menopause_status'),\n",
    "    ('birth_days_to', 'birth_days_to'),\n",
    "    ('vital_status', 'vital_status'),\n",
    "    ('tumor_status', 'tumor_status'),\n",
    "    ('last_contact_days_to', 'last_contact_days_to'),\n",
    "    ('death_days_to', 'death_days_to'),\n",
    "    ('cause_of_death', 'cause_of_death'),\n",
    "    ('new_tumor_event_type', 'new_tumor_event_type'),\n",
    "    ('new_tumor_event_site', 'new_tumor_event_site'),\n",
    "    ('new_tumor_event_site_other', 'new_tumor_event_site_other'),\n",
    "    ('new_tumor_event_dx_days_to', 'days_recurrence_free'),\n",
    "    ('treatment_outcome_first_course', 'treatment_outcome_first_course'),\n",
    "    ('margin_status', 'margin_status'),\n",
    "    ('residual_tumor', 'residual_tumor'),\n",
    "    ('_EVENT', 'event_status'),\n",
    "    ('_TIME_TO_EVENT', 'event_days'),\n",
    "    ('OS', 'dead'),\n",
    "    ('OS.time', 'days_survived'),\n",
    "    ('DSS', 'disease_specific_survival_status'),\n",
    "    ('DSS.time', 'disease_specific_survival_days'),\n",
    "    ('DFI', 'disease_free_interval_status'),\n",
    "    ('DFI.time', 'disease_free_interval_days'),\n",
    "    ('PFI', 'progression_free_interval_status'),\n",
    "    ('PFI.time', 'progression_free_interval_days'),\n",
    "    ('Redaction', 'Redaction')\n",
    "])\n",
    "\n",
    "clinmat_df = (\n",
    "    pandas.read_table(path)\n",
    "    .rename(columns=renamer)\n",
    "    .merge(disease_df, how='left')\n",
    "    [list(renamer.values())]\n",
    "    .set_index('sample_id', drop=False)\n",
    "    .sort_values('sample_id')\n",
    ")\n",
    "\n",
    "# Fix capitalization of gender and race\n",
    "clinmat_df.gender = clinmat_df.gender.str.title()\n",
    "clinmat_df.race = clinmat_df.race.str.title()\n",
    "\n",
    "# Extract sample-type with the code dictionary\n",
    "clinmat_df = clinmat_df.assign(sample_type = clinmat_df.sample_id.str[-2:])\n",
    "clinmat_df.sample_type = clinmat_df.sample_type.replace(sampletype_codes_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sample_id</th>\n",
       "      <th>patient_id</th>\n",
       "      <th>acronym</th>\n",
       "      <th>disease</th>\n",
       "      <th>age_diagnosed</th>\n",
       "      <th>gender</th>\n",
       "      <th>race</th>\n",
       "      <th>ajcc_stage</th>\n",
       "      <th>clinical_stage</th>\n",
       "      <th>histological_type</th>\n",
       "      <th>...</th>\n",
       "      <th>dead</th>\n",
       "      <th>days_survived</th>\n",
       "      <th>disease_specific_survival_status</th>\n",
       "      <th>disease_specific_survival_days</th>\n",
       "      <th>disease_free_interval_status</th>\n",
       "      <th>disease_free_interval_days</th>\n",
       "      <th>progression_free_interval_status</th>\n",
       "      <th>progression_free_interval_days</th>\n",
       "      <th>Redaction</th>\n",
       "      <th>sample_type</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0001-01</th>\n",
       "      <td>TCGA-02-0001-01</td>\n",
       "      <td>TCGA-02-0001</td>\n",
       "      <td>GBM</td>\n",
       "      <td>glioblastoma multiforme</td>\n",
       "      <td>44.0</td>\n",
       "      <td>Female</td>\n",
       "      <td>White</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Untreated primary (de novo) GBM</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>358.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>358.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>137.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Primary Solid Tumor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0003-01</th>\n",
       "      <td>TCGA-02-0003-01</td>\n",
       "      <td>TCGA-02-0003</td>\n",
       "      <td>GBM</td>\n",
       "      <td>glioblastoma multiforme</td>\n",
       "      <td>50.0</td>\n",
       "      <td>Male</td>\n",
       "      <td>White</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Untreated primary (de novo) GBM</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Primary Solid Tumor</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 38 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                       sample_id    patient_id acronym  \\\n",
       "sample_id                                                \n",
       "TCGA-02-0001-01  TCGA-02-0001-01  TCGA-02-0001     GBM   \n",
       "TCGA-02-0003-01  TCGA-02-0003-01  TCGA-02-0003     GBM   \n",
       "\n",
       "                                 disease  age_diagnosed  gender   race  \\\n",
       "sample_id                                                                \n",
       "TCGA-02-0001-01  glioblastoma multiforme           44.0  Female  White   \n",
       "TCGA-02-0003-01  glioblastoma multiforme           50.0    Male  White   \n",
       "\n",
       "                ajcc_stage clinical_stage                histological_type  \\\n",
       "sample_id                                                                    \n",
       "TCGA-02-0001-01        NaN            NaN  Untreated primary (de novo) GBM   \n",
       "TCGA-02-0003-01        NaN            NaN  Untreated primary (de novo) GBM   \n",
       "\n",
       "                        ...          dead  days_survived  \\\n",
       "sample_id               ...                                \n",
       "TCGA-02-0001-01         ...           1.0          358.0   \n",
       "TCGA-02-0003-01         ...           1.0          144.0   \n",
       "\n",
       "                disease_specific_survival_status  \\\n",
       "sample_id                                          \n",
       "TCGA-02-0001-01                              1.0   \n",
       "TCGA-02-0003-01                              1.0   \n",
       "\n",
       "                 disease_specific_survival_days disease_free_interval_status  \\\n",
       "sample_id                                                                      \n",
       "TCGA-02-0001-01                           358.0                          NaN   \n",
       "TCGA-02-0003-01                           144.0                          NaN   \n",
       "\n",
       "                disease_free_interval_days  progression_free_interval_status  \\\n",
       "sample_id                                                                      \n",
       "TCGA-02-0001-01                        NaN                               1.0   \n",
       "TCGA-02-0003-01                        NaN                               1.0   \n",
       "\n",
       "                 progression_free_interval_days Redaction          sample_type  \n",
       "sample_id                                                                       \n",
       "TCGA-02-0001-01                           137.0       NaN  Primary Solid Tumor  \n",
       "TCGA-02-0003-01                            40.0       NaN  Primary Solid Tumor  \n",
       "\n",
       "[2 rows x 38 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clinmat_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Primary Solid Tumor                                10517\n",
       "Solid Tissue Normal                                 1413\n",
       "Metastatic                                           395\n",
       "Primary Blood Derived Cancer - Peripheral Blood      200\n",
       "Recurrent Solid Tumor                                 55\n",
       "Additional - New Primary                              10\n",
       "Additional Metastatic                                  1\n",
       "Name: sample_type, dtype: int64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What is the distribution of cancer-types\n",
    "clinmat_df.sample_type.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save unfiltered dataset to a TSV\n",
    "path = os.path.join('data', 'complete', 'samples.tsv')\n",
    "clinmat_df.to_csv(path, sep='\\t', float_format='%.0f', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove \"Redacted\" tumors\n",
    "# These patients either withdrew consent or had genome data mismatch errors\n",
    "clinmat_df = clinmat_df.query('Redaction != \"Redacted\"').drop(\"Redaction\", axis=1)\n",
    "\n",
    "# Keep only these sample types\n",
    "# See https://github.com/cognoma/cancer-data/issues/46\n",
    "sample_types = {\n",
    "    'Primary Solid Tumor',\n",
    "    'Recurrent Solid Tumor',\n",
    "    'Primary Blood Derived Cancer - Peripheral Blood',\n",
    "    'Metastatic',\n",
    "}\n",
    "\n",
    "clinmat_df = (\n",
    "    clinmat_df\n",
    "    .query(\"sample_type in @sample_types\")\n",
    "    # filters duplicate samples per patient (sample_ids are sorted,\n",
    "    # so two-digit suffix places a metastatic tumor after a primary\n",
    "    # solid tumor for a given a patient).\n",
    "    .drop_duplicates('patient_id', keep='first')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11026"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check that no samples are duplicated\n",
    "assert not clinmat_df.duplicated('sample_id', keep=False).any()\n",
    "\n",
    "# Check that all diseases in clinmat_df are in disease_df\n",
    "assert not set(clinmat_df.acronym) - set(disease_df.acronym)\n",
    "\n",
    "len(clinmat_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sample_id</th>\n",
       "      <th>patient_id</th>\n",
       "      <th>acronym</th>\n",
       "      <th>disease</th>\n",
       "      <th>age_diagnosed</th>\n",
       "      <th>gender</th>\n",
       "      <th>race</th>\n",
       "      <th>ajcc_stage</th>\n",
       "      <th>clinical_stage</th>\n",
       "      <th>histological_type</th>\n",
       "      <th>...</th>\n",
       "      <th>event_days</th>\n",
       "      <th>dead</th>\n",
       "      <th>days_survived</th>\n",
       "      <th>disease_specific_survival_status</th>\n",
       "      <th>disease_specific_survival_days</th>\n",
       "      <th>disease_free_interval_status</th>\n",
       "      <th>disease_free_interval_days</th>\n",
       "      <th>progression_free_interval_status</th>\n",
       "      <th>progression_free_interval_days</th>\n",
       "      <th>sample_type</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0001-01</th>\n",
       "      <td>TCGA-02-0001-01</td>\n",
       "      <td>TCGA-02-0001</td>\n",
       "      <td>GBM</td>\n",
       "      <td>glioblastoma multiforme</td>\n",
       "      <td>44.0</td>\n",
       "      <td>Female</td>\n",
       "      <td>White</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Untreated primary (de novo) GBM</td>\n",
       "      <td>...</td>\n",
       "      <td>358.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>358.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>358.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>137.0</td>\n",
       "      <td>Primary Solid Tumor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0003-01</th>\n",
       "      <td>TCGA-02-0003-01</td>\n",
       "      <td>TCGA-02-0003</td>\n",
       "      <td>GBM</td>\n",
       "      <td>glioblastoma multiforme</td>\n",
       "      <td>50.0</td>\n",
       "      <td>Male</td>\n",
       "      <td>White</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Untreated primary (de novo) GBM</td>\n",
       "      <td>...</td>\n",
       "      <td>144.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>Primary Solid Tumor</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 37 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                       sample_id    patient_id acronym  \\\n",
       "sample_id                                                \n",
       "TCGA-02-0001-01  TCGA-02-0001-01  TCGA-02-0001     GBM   \n",
       "TCGA-02-0003-01  TCGA-02-0003-01  TCGA-02-0003     GBM   \n",
       "\n",
       "                                 disease  age_diagnosed  gender   race  \\\n",
       "sample_id                                                                \n",
       "TCGA-02-0001-01  glioblastoma multiforme           44.0  Female  White   \n",
       "TCGA-02-0003-01  glioblastoma multiforme           50.0    Male  White   \n",
       "\n",
       "                ajcc_stage clinical_stage                histological_type  \\\n",
       "sample_id                                                                    \n",
       "TCGA-02-0001-01        NaN            NaN  Untreated primary (de novo) GBM   \n",
       "TCGA-02-0003-01        NaN            NaN  Untreated primary (de novo) GBM   \n",
       "\n",
       "                        ...          event_days  dead days_survived  \\\n",
       "sample_id               ...                                           \n",
       "TCGA-02-0001-01         ...               358.0   1.0         358.0   \n",
       "TCGA-02-0003-01         ...               144.0   1.0         144.0   \n",
       "\n",
       "                 disease_specific_survival_status  \\\n",
       "sample_id                                           \n",
       "TCGA-02-0001-01                               1.0   \n",
       "TCGA-02-0003-01                               1.0   \n",
       "\n",
       "                disease_specific_survival_days disease_free_interval_status  \\\n",
       "sample_id                                                                     \n",
       "TCGA-02-0001-01                          358.0                          NaN   \n",
       "TCGA-02-0003-01                          144.0                          NaN   \n",
       "\n",
       "                 disease_free_interval_days  progression_free_interval_status  \\\n",
       "sample_id                                                                       \n",
       "TCGA-02-0001-01                         NaN                               1.0   \n",
       "TCGA-02-0003-01                         NaN                               1.0   \n",
       "\n",
       "                progression_free_interval_days          sample_type  \n",
       "sample_id                                                            \n",
       "TCGA-02-0001-01                          137.0  Primary Solid Tumor  \n",
       "TCGA-02-0003-01                           40.0  Primary Solid Tumor  \n",
       "\n",
       "[2 rows x 37 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clinmat_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read mutation data\n",
    "\n",
    "This file contains mutation data (which mutations each sample contains). See the [online documentation](https://xenabrowser.net/datapages/?dataset=mc3.v0.2.8.PUBLIC.xena&host=https%3A%2F%2Fpancanatlas.xenahubs.net) for the `MC3` resource. For more information about mutation calling refer to [Ellrott et al. 2018](https://doi.org/10.1016/j.cels.2018.03.002 \"Scalable Open Science Approach for Mutation Calling of Tumor Exomes Using Multiple Genomic Pipelines\")\n",
    "\n",
    "Note that duplicate mutation rows, which [occur](https://groups.google.com/d/msg/ucsc-cancer-genomics-browser/eg6nJOFSefw/Z0BM6pU9BAAJ \"Message on the Xena Browser Google Group\") for samples that were sequenced multiple times, are filtered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>symbol</th>\n",
       "      <th>chromosome</th>\n",
       "      <th>entrez_gene_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>(FM-3)</td>\n",
       "      <td>chr2</td>\n",
       "      <td>10316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>(IV)-44</td>\n",
       "      <td>chr14</td>\n",
       "      <td>28337</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    symbol chromosome  entrez_gene_id\n",
       "0   (FM-3)       chr2           10316\n",
       "1  (IV)-44      chr14           28337"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load chromosome-symbol to entrez_gene_id mapping\n",
    "path = os.path.join('download', 'genes', 'chromosome-symbol-mapper.tsv')\n",
    "chr_sym_map_df = pandas.read_table(path)\n",
    "chr_sym_map_df.chromosome = 'chr' + chr_sym_map_df.chromosome\n",
    "chr_sym_map_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sample_id</th>\n",
       "      <th>chr</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>reference</th>\n",
       "      <th>alt</th>\n",
       "      <th>effect</th>\n",
       "      <th>Amino_Acid_Change</th>\n",
       "      <th>DNA_VAF</th>\n",
       "      <th>SIFT</th>\n",
       "      <th>PolyPhen</th>\n",
       "      <th>entrez_gene_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TCGA-02-0003-01</td>\n",
       "      <td>10</td>\n",
       "      <td>123810032</td>\n",
       "      <td>123810032</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>p.T38M</td>\n",
       "      <td>0.88</td>\n",
       "      <td>NaN</td>\n",
       "      <td>benign(0.335)</td>\n",
       "      <td>10579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TCGA-05-4427-01</td>\n",
       "      <td>10</td>\n",
       "      <td>123846949</td>\n",
       "      <td>123846949</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>p.E1645V</td>\n",
       "      <td>0.49</td>\n",
       "      <td>NaN</td>\n",
       "      <td>possibly_damaging(0.489)</td>\n",
       "      <td>10579</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         sample_id chr      start        end reference alt             effect  \\\n",
       "0  TCGA-02-0003-01  10  123810032  123810032         C   T  Missense_Mutation   \n",
       "1  TCGA-05-4427-01  10  123846949  123846949         A   T  Missense_Mutation   \n",
       "\n",
       "  Amino_Acid_Change  DNA_VAF SIFT                  PolyPhen  entrez_gene_id  \n",
       "0            p.T38M     0.88  NaN             benign(0.335)           10579  \n",
       "1          p.E1645V     0.49  NaN  possibly_damaging(0.489)           10579  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = os.path.join('download', 'mc3.v0.2.8.PUBLIC.xena.tsv.gz')\n",
    "snp_mutation_df = (\n",
    "    pandas.read_table(path)\n",
    "    .rename(columns={'sample': 'sample_id', 'gene': 'symbol'})\n",
    "    .merge(chr_sym_map_df, on='symbol')\n",
    "    .drop(['chromosome', 'symbol'], axis='columns')\n",
    "    .drop_duplicates()\n",
    ")\n",
    "snp_mutation_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9104"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of samples with at least one mutation\n",
    "samples_with_mutation_calls = sorted(set(snp_mutation_df.sample_id))\n",
    "len(samples_with_mutation_calls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20393"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of genes with at least one mutation\n",
    "snp_mutation_df.entrez_gene_id.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>effect</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Missense_Mutation</td>\n",
       "      <td>1606302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Silent</td>\n",
       "      <td>654079</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3'UTR</td>\n",
       "      <td>251556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Nonsense_Mutation</td>\n",
       "      <td>126716</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Intron</td>\n",
       "      <td>88756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Frame_Shift_Del</td>\n",
       "      <td>77710</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>5'UTR</td>\n",
       "      <td>69267</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Splice_Site</td>\n",
       "      <td>40180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>RNA</td>\n",
       "      <td>34116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Frame_Shift_Ins</td>\n",
       "      <td>22668</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>3'Flank</td>\n",
       "      <td>17240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5'Flank</td>\n",
       "      <td>10704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>In_Frame_Del</td>\n",
       "      <td>8939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Translation_Start_Site</td>\n",
       "      <td>2235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Nonstop_Mutation</td>\n",
       "      <td>1677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>In_Frame_Ins</td>\n",
       "      <td>789</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>large deletion</td>\n",
       "      <td>159</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     index   effect\n",
       "0        Missense_Mutation  1606302\n",
       "1                   Silent   654079\n",
       "2                    3'UTR   251556\n",
       "3        Nonsense_Mutation   126716\n",
       "4                   Intron    88756\n",
       "5          Frame_Shift_Del    77710\n",
       "6                    5'UTR    69267\n",
       "7              Splice_Site    40180\n",
       "8                      RNA    34116\n",
       "9          Frame_Shift_Ins    22668\n",
       "10                 3'Flank    17240\n",
       "11                 5'Flank    10704\n",
       "12            In_Frame_Del     8939\n",
       "13  Translation_Start_Site     2235\n",
       "14        Nonstop_Mutation     1677\n",
       "15            In_Frame_Ins      789\n",
       "16          large deletion      159"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Mutations counts by type\n",
    "snp_mutation_df.effect.value_counts().reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert SNP mutations to gene mutations\n",
    "\n",
    "The next cell specifies which mutations to preserve as gene-affecting, which were chosen according to the red & blue [mutation effects in Xena](http://xena.ucsc.edu/how-we-characterize-mutations/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "mutations = {\n",
    "    'Frame_Shift_Del',\n",
    "    'Frame_Shift_Ins',\n",
    "    'In_Frame_Del',\n",
    "    'In_Frame_Ins',\n",
    "    'Missense_Mutation',\n",
    "    'Nonsense_Mutation',\n",
    "    'Nonstop_Mutation',\n",
    "    'RNA',\n",
    "    'Splice_Site',\n",
    "    'Translation_Start_Site',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{\"3'Flank\", \"3'UTR\", \"5'Flank\", \"5'UTR\", 'Intron', 'Silent', 'large deletion'}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Mutations effects that were observed but nut included\n",
    "set(snp_mutation_df.effect.unique()) - mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sample_id</th>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>TCGA-02-0003-01</td>\n",
       "      <td>1665</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>TCGA-02-0003-01</td>\n",
       "      <td>1956</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         sample_id  entrez_gene_id  count\n",
       "0  TCGA-02-0003-01            1665      1\n",
       "1  TCGA-02-0003-01            1956      1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_mutation_df = (snp_mutation_df\n",
    "    .query(\"effect in @mutations\")\n",
    "    .groupby(['sample_id', 'entrez_gene_id'])\n",
    "    .apply(len)\n",
    "    .reset_index()\n",
    "    .rename(columns={0: 'count'})\n",
    ")\n",
    "\n",
    "gene_mutation_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9104, 20224)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a sample (rows) by gene (columns) matrix of mutation status\n",
    "gene_mutation_mat_df = (gene_mutation_df\n",
    "    .pivot_table(index='sample_id',\n",
    "                 columns='entrez_gene_id',\n",
    "                 values='count',\n",
    "                 fill_value=0)\n",
    "    .reindex(samples_with_mutation_calls, fill_value=0)\n",
    "    .astype(bool).astype(int)\n",
    ")\n",
    "gene_mutation_mat_df.columns = gene_mutation_mat_df.columns.astype(str)\n",
    "gene_mutation_mat_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.90% sample-gene pairs are mutated'"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'{:.2%} sample-gene pairs are mutated'.format(\n",
    "    gene_mutation_mat_df.stack().mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save complete mutation matrix\n",
    "path = os.path.join('data', 'complete', 'mutation-matrix.tsv.bz2')\n",
    "gene_mutation_mat_df.to_csv(path, sep='\\t', compression='bz2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read gene expression data\n",
    "\n",
    "This file contains gene expression data. See the [online documentation](https://xenabrowser.net/datapages/?dataset=EB%2B%2BAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.xena&host=https%3A%2F%2Fpancanatlas.xenahubs.net) for processing information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the gene × sample dataset\n",
    "path = os.path.join('download', 'EB++AdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.xena.tsv.gz')\n",
    "expr_df = pandas.read_table(path, index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'100130426',\n",
       " '100133144',\n",
       " '100134869',\n",
       " '10357',\n",
       " '10431',\n",
       " '136542',\n",
       " '155060',\n",
       " '26823',\n",
       " '280660',\n",
       " '317712',\n",
       " '340602',\n",
       " '388795',\n",
       " '390284',\n",
       " '391343',\n",
       " '391714',\n",
       " '404770',\n",
       " '441362',\n",
       " '442388',\n",
       " '553137',\n",
       " '57714',\n",
       " '645851',\n",
       " '652919',\n",
       " '653553',\n",
       " '728045',\n",
       " '728603',\n",
       " '728788',\n",
       " '729884',\n",
       " '8225',\n",
       " '90288'}"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Retrieve symbol to gene mapping for gene expression data\n",
    "path = os.path.join('mapping', 'HiSeqV2-genes', 'HiSeqV2-gene-map.tsv')\n",
    "gene_map_df = pandas.read_table(path, dtype='str')\n",
    "symbol_to_entrez = dict(zip(gene_map_df.symbol, gene_map_df.entrez_gene_id))\n",
    "\n",
    "# Check for unmapped symbols\n",
    "unmapped_symbols = set(expr_df.index) - set(symbol_to_entrez)\n",
    "unmapped_symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11069, 16261)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Process the dataset\n",
    "expr_df = (expr_df\n",
    "    # Drop genes with NA measurements\n",
    "    # (See https://github.com/cognoma/cancer-data/issues/40)\n",
    "    .dropna(axis='rows')\n",
    "    # Convert gene symbols to entrez gene ids\n",
    "    .rename(index=symbol_to_entrez)\n",
    "    .rename(index=old_to_new_entrez)\n",
    "    # Average expression for rows with the same entrez_gene_id\n",
    "    .groupby(level=0).mean()\n",
    "    # Transpose so the data is sample × gene\n",
    "    .transpose()\n",
    "    # Sort rows and columns\n",
    "    .sort_index(axis='rows')\n",
    "    .sort_index(axis='columns')\n",
    ")\n",
    "\n",
    "expr_df.index.rename('sample_id', inplace=True)\n",
    "\n",
    "# Filter for valid Entrez Genes (This will filter the unmapped symbols)\n",
    "expr_df = expr_df.loc[:, expr_df.columns.isin(gene_df.entrez_gene_id)]\n",
    "\n",
    "expr_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10193"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of patients represented in the expression dataset\n",
    "clinmat_df.query(\"sample_id in @expr_df.index\").patient_id.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sample</th>\n",
       "      <th>1</th>\n",
       "      <th>100</th>\n",
       "      <th>1000</th>\n",
       "      <th>10000</th>\n",
       "      <th>100009676</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0047-01</th>\n",
       "      <td>6.98</td>\n",
       "      <td>7.10</td>\n",
       "      <td>11.17</td>\n",
       "      <td>10.34</td>\n",
       "      <td>6.18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-0055-01</th>\n",
       "      <td>8.62</td>\n",
       "      <td>7.80</td>\n",
       "      <td>10.83</td>\n",
       "      <td>9.82</td>\n",
       "      <td>6.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-2483-01</th>\n",
       "      <td>8.09</td>\n",
       "      <td>8.00</td>\n",
       "      <td>11.50</td>\n",
       "      <td>10.37</td>\n",
       "      <td>6.12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-2485-01</th>\n",
       "      <td>6.41</td>\n",
       "      <td>7.02</td>\n",
       "      <td>12.77</td>\n",
       "      <td>13.31</td>\n",
       "      <td>5.18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCGA-02-2486-01</th>\n",
       "      <td>6.77</td>\n",
       "      <td>7.69</td>\n",
       "      <td>11.14</td>\n",
       "      <td>9.77</td>\n",
       "      <td>6.24</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "sample              1   100   1000  10000  100009676\n",
       "sample_id                                           \n",
       "TCGA-02-0047-01  6.98  7.10  11.17  10.34       6.18\n",
       "TCGA-02-0055-01  8.62  7.80  10.83   9.82       6.60\n",
       "TCGA-02-2483-01  8.09  8.00  11.50  10.37       6.12\n",
       "TCGA-02-2485-01  6.41  7.02  12.77  13.31       5.18\n",
       "TCGA-02-2486-01  6.77  7.69  11.14   9.77       6.24"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Peak at the data matrix\n",
    "expr_df.iloc[:5, :5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save complete expression matrix\n",
    "path = os.path.join('data', 'complete', 'expression-matrix.tsv.bz2')\n",
    "expr_df.to_csv(path, sep='\\t', float_format='%.3g', compression='bz2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrate expression and mutation data\n",
    "\n",
    "Find samples with both mutation and expression data.\n",
    "\n",
    "We assume that if a sample was not in the `MC3` data, it was not assayed for mutation ([more info](https://github.com/cognoma/cancer-data/issues/43#issuecomment-380957274))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8756"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_ids = sorted(clinmat_df.index & gene_mutation_mat_df.index & expr_df.index)\n",
    "len(sample_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter expression (x) and mutation (y) matrices for common samples\n",
    "sample_df = clinmat_df.loc[sample_ids, :]\n",
    "x_df = expr_df.loc[sample_ids, :]\n",
    "y_df = gene_mutation_mat_df.loc[sample_ids, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a columnn for mutations per sample\n",
    "sample_df['n_mutations'] = y_df.sum(axis='columns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th>symbol</th>\n",
       "      <th>description</th>\n",
       "      <th>gene_type</th>\n",
       "      <th>mean_expression</th>\n",
       "      <th>stdev_expression</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>A1BG</td>\n",
       "      <td>alpha-1-B glycoprotein</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>6.591484</td>\n",
       "      <td>2.281684</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>100</td>\n",
       "      <td>ADA</td>\n",
       "      <td>adenosine deaminase</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>7.474607</td>\n",
       "      <td>1.525839</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       entrez_gene_id symbol             description       gene_type  \\\n",
       "sample                                                                 \n",
       "1                   1   A1BG  alpha-1-B glycoprotein  protein-coding   \n",
       "100               100    ADA     adenosine deaminase  protein-coding   \n",
       "\n",
       "        mean_expression  stdev_expression  \n",
       "sample                                     \n",
       "1              6.591484          2.281684  \n",
       "100            7.474607          1.525839  "
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_gene_df = (\n",
    "    gene_df.loc[x_df.columns, :]\n",
    "    .assign(mean_expression=x_df.mean(axis='rows'))\n",
    "    .assign(stdev_expression=x_df.std(axis='rows'))\n",
    ")\n",
    "path = os.path.join('data', 'expression-genes.tsv')\n",
    "x_gene_df.to_csv(path, sep='\\t', index=False, float_format='%.4g')\n",
    "x_gene_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th>symbol</th>\n",
       "      <th>description</th>\n",
       "      <th>gene_type</th>\n",
       "      <th>n_mutations</th>\n",
       "      <th>mutation_freq</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>entrez_gene_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>A1BG</td>\n",
       "      <td>alpha-1-B glycoprotein</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>57</td>\n",
       "      <td>0.006510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>A2M</td>\n",
       "      <td>alpha-2-macroglobulin</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>226</td>\n",
       "      <td>0.025811</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               entrez_gene_id symbol             description       gene_type  \\\n",
       "entrez_gene_id                                                                 \n",
       "1                           1   A1BG  alpha-1-B glycoprotein  protein-coding   \n",
       "2                           2    A2M   alpha-2-macroglobulin  protein-coding   \n",
       "\n",
       "                n_mutations  mutation_freq  \n",
       "entrez_gene_id                              \n",
       "1                        57       0.006510  \n",
       "2                       226       0.025811  "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_gene_df = (\n",
    "    gene_df.loc[y_df.columns, :]\n",
    "    .assign(n_mutations=y_df.sum(axis='rows'))\n",
    "    .assign(mutation_freq=y_df.mean(axis='rows'))\n",
    ")\n",
    "path = os.path.join('data', 'mutation-genes.tsv')\n",
    "y_gene_df.to_csv(path, sep='\\t', index=False, float_format='%.4g')\n",
    "y_gene_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cancer type (disease) stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>acronym</th>\n",
       "      <th>disease</th>\n",
       "      <th>n_samples</th>\n",
       "      <th>n_clinical_samples</th>\n",
       "      <th>n_expression_samples</th>\n",
       "      <th>n_mutation_samples</th>\n",
       "      <th>median_mutations</th>\n",
       "      <th>mean_mutations</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ACC</td>\n",
       "      <td>adrenocortical cancer</td>\n",
       "      <td>79</td>\n",
       "      <td>92</td>\n",
       "      <td>79</td>\n",
       "      <td>92</td>\n",
       "      <td>27.0</td>\n",
       "      <td>68.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BLCA</td>\n",
       "      <td>bladder urothelial carcinoma</td>\n",
       "      <td>403</td>\n",
       "      <td>409</td>\n",
       "      <td>405</td>\n",
       "      <td>407</td>\n",
       "      <td>164.0</td>\n",
       "      <td>228.168734</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  acronym                       disease  n_samples  n_clinical_samples  \\\n",
       "0     ACC         adrenocortical cancer         79                  92   \n",
       "1    BLCA  bladder urothelial carcinoma        403                 409   \n",
       "\n",
       "   n_expression_samples  n_mutation_samples  median_mutations  mean_mutations  \n",
       "0                    79                  92              27.0       68.000000  \n",
       "1                   405                 407             164.0      228.168734  "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_to_acronym = dict(zip(clinmat_df.sample_id, clinmat_df.acronym))\n",
    "\n",
    "def get_cancer_count_column(sample_ids):\n",
    "    \"\"\"\n",
    "    sample_ids is a pandas.Series\n",
    "    \"\"\"\n",
    "    sample_ids = pandas.Series(sample_ids)\n",
    "    acronyms = sample_ids.map(sample_to_acronym)\n",
    "    counter = collections.Counter(acronyms)\n",
    "    counts = disease_df.acronym.map(counter)\n",
    "    return counts.fillna(0).astype(int)\n",
    "\n",
    "# Compute nubmer of samples per disease (cancer type)\n",
    "disease_df['n_samples'] = get_cancer_count_column(sample_df.sample_id)\n",
    "disease_df['n_clinical_samples'] = get_cancer_count_column(clinmat_df.sample_id)\n",
    "disease_df['n_expression_samples'] = get_cancer_count_column(expr_df.index)\n",
    "disease_df['n_mutation_samples'] = get_cancer_count_column(gene_mutation_mat_df.index)\n",
    "\n",
    "# Compute n_mutation summaries for samples in the aligned set\n",
    "acronyms = list(pandas.Series(y_df.index).map(sample_to_acronym))\n",
    "groups = y_df.sum(axis='columns').groupby(acronyms)\n",
    "disease_df['median_mutations'] = disease_df.acronym.map(dict(groups.median()))\n",
    "disease_df['mean_mutations'] = disease_df.acronym.map(dict(groups.mean()))\n",
    "\n",
    "# Export to TSV\n",
    "path = os.path.join('data', 'diseases.tsv')\n",
    "disease_df.to_csv(path, sep='\\t', float_format='%.1f', index=False)\n",
    "\n",
    "disease_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export matrices to TSVs\n",
    "\n",
    "Matrices are saved as sample × gene TSVs. Subsetted matrices are also exported to allow users to quickly explore small portions of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = os.path.join('data', 'samples.tsv')\n",
    "sample_df.to_csv(path, sep='\\t', float_format='%.0f', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subset_df(df, nrows=None, ncols=None, row_seed=0, col_seed=0):\n",
    "    \"\"\"Randomly subset a dataframe, preserving row and column order.\"\"\"\n",
    "    if nrows is None:\n",
    "        nrows = len(df)\n",
    "    if ncols is None:\n",
    "        ncols = len(df.columns)\n",
    "    return (df\n",
    "        .sample(n=nrows, random_state=row_seed, axis='rows')\n",
    "        .sample(n=ncols, random_state=col_seed, axis='columns')\n",
    "        .sort_index(axis='rows')\n",
    "        .sort_index(axis='columns')\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsv_args = {'sep': '\\t', 'float_format': '%.3g'}\n",
    "\n",
    "for df, name in (x_df, 'expression-matrix'), (y_df, 'mutation-matrix'):\n",
    "\n",
    "    # Save full dataset\n",
    "    path = os.path.join('data', name + '.tsv.bz2')\n",
    "    df.to_csv(path, **tsv_args, compression='bz2')\n",
    "    \n",
    "    # Save subsetted datasets\n",
    "    for sample, nrows, ncols in ('small', 50, 15), ('all-samples', None, 15), ('all-genes', 50, None):\n",
    "        path = os.path.join('data', 'subset', '{}-{}.tsv'.format(name, sample))\n",
    "        subset_df(df, nrows=nrows, ncols=ncols).to_csv(path, **tsv_args)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:cognoma-cancer-data]",
   "language": "python",
   "name": "conda-env-cognoma-cancer-data-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}