survival_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TGCA Survival Analysis of Pancreatic Cancer Patients\n",
    "\n",
    "### Modified from similar workflow described on biostars.org: https://www.biostars.org/p/153013/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "from collections import Counter\n",
    "import pandas as pd\n",
    "import rpy2.rinterface\n",
    "import numpy as np\n",
    "import glob\n",
    "import os\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download raw counts from Xena, place tsv in data folder\n",
    "# https://gdc.xenahubs.net/download/TCGA-PAAD/Xena_Matrices/TCGA-PAAD.htseq_counts.tsv.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "HTSeq_raw_df = pd.read_csv(\"data/rna/raw_counts/TCGA-PAAD.htseq_counts.tsv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These counts are actually partially processed from raw HTSeq counts (they have been incremented by one then log transformed)\n",
    "# Reversing this process will give the actual raw counts\n",
    "\n",
    "HTSeq_raw_counts_df = pd.DataFrame()\n",
    "HTSeq_raw_counts_df['Ensembl_ID'] = HTSeq_raw_df['Ensembl_ID']\n",
    "for column in HTSeq_raw_df.columns.values[1:]:\n",
    "    series = HTSeq_raw_df[column]\n",
    "    series_antilog = (2**series) - 1\n",
    "    HTSeq_raw_counts_df[column] = series_antilog\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Builds a dictionary of what samples correspond to normal tissue and what samples correspond to tumor tissue\n",
    "# More can be read about how sample barcoding works in the TCGA here: https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/\n",
    "\n",
    "origin_dict = {}\n",
    "for column in list(HTSeq_raw_counts_df)[1:]:\n",
    "    if column[13] == '0':\n",
    "        origin_dict[column] = 'tumor'\n",
    "    elif column[13] == '1':\n",
    "        origin_dict[column] = 'normal'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({'tumor': 178, 'normal': 4})"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# There should be 178 samples corresponding to tumors and 4 samples corresponding to normal tissue\n",
    "\n",
    "Counter(origin_dict.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For this analysis, only tumor samples are used\n",
    "# HTSeq_raw_counts_df is subset to just get tumor samples\n",
    "# Normal samples are not used in this analysis\n",
    "\n",
    "slice_list = ['Ensembl_ID']\n",
    "for column, value in origin_dict.items():\n",
    "    if value == 'tumor':\n",
    "        slice_list.append(column)\n",
    "        \n",
    "HTSeq_raw_counts_tumors_df = HTSeq_raw_counts_df[slice_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lowly expressed genes (genes that have counts of 0 in 50% of samples) are removed\n",
    "\n",
    "zeroes = (HTSeq_raw_counts_tumors_df == 0).astype(int)\n",
    "zeroes_count = zeroes.sum(axis=1).to_frame()\n",
    "length = zeroes.shape[1]\n",
    "zeroes_count = zeroes_count.drop(zeroes_count[zeroes_count[0] < length/2].index)\n",
    "drop_index_list = zeroes_count.index.to_list()\n",
    "HTSeq_raw_counts_tumors_no_zeroes_df = HTSeq_raw_counts_tumors_df.drop(labels=drop_index_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Counts are floats currently, but need to be convereted to integers of DESeq2 will complain\n",
    "\n",
    "for column in HTSeq_raw_counts_tumors_no_zeroes_df.columns.values[1:]:\n",
    "    HTSeq_raw_counts_tumors_no_zeroes_df[column] = HTSeq_raw_counts_tumors_no_zeroes_df[column].astype(int)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Any files that are currently in the tumor_samples folder should be deleted\n",
    "\n",
    "tumor_sample_dir = 'data/rna/tumor_samples/'\n",
    "files = glob.glob(tumor_sample_dir + '*')\n",
    "for f in files:\n",
    "    os.remove(f)\n",
    "    \n",
    "# A tsv corresponing to the raw counts for each tumor is saved in the tumor_samples folder\n",
    "# This is needed to pass off the data to R\n",
    "\n",
    "for column in HTSeq_raw_counts_tumors_no_zeroes_df.columns.values[1:]:\n",
    "    file_name = 'data/rna/tumor_samples/' + column + '.tsv'\n",
    "    HTSeq_raw_counts_tumors_no_zeroes_df[['Ensembl_ID', column]].to_csv(file_name, sep=\" \", header=False, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: S4Vectors\n",
      "\n",
      "R[write to console]: Loading required package: stats4\n",
      "\n",
      "R[write to console]: Loading required package: BiocGenerics\n",
      "\n",
      "R[write to console]: Loading required package: parallel\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:parallel’:\n",
      "\n",
      "    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
      "    clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
      "    parLapplyLB, parRapply, parSapply, parSapplyLB\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, sd, var, xtabs\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, append, as.data.frame, basename, cbind, colnames,\n",
      "    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,\n",
      "    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,\n",
      "    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,\n",
      "    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,\n",
      "    union, unique, unsplit, which, which.max, which.min\n",
      "\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘S4Vectors’\n",
      "\n",
      "\n",
      "R[write to console]: The following object is masked from ‘package:base’:\n",
      "\n",
      "    expand.grid\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: IRanges\n",
      "\n",
      "R[write to console]: Loading required package: GenomicRanges\n",
      "\n",
      "R[write to console]: Loading required package: GenomeInfoDb\n",
      "\n",
      "R[write to console]: Loading required package: SummarizedExperiment\n",
      "\n",
      "R[write to console]: Loading required package: Biobase\n",
      "\n",
      "R[write to console]: Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: DelayedArray\n",
      "\n",
      "R[write to console]: Loading required package: matrixStats\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘matrixStats’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:Biobase’:\n",
      "\n",
      "    anyMissing, rowMedians\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: BiocParallel\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘DelayedArray’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:matrixStats’:\n",
      "\n",
      "    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    aperm, apply, rowsum\n",
      "\n",
      "\n",
      "R[write to console]: Registered S3 methods overwritten by 'ggplot2':\n",
      "  method         from \n",
      "  [.quosures     rlang\n",
      "  c.quosures     rlang\n",
      "  print.quosures rlang\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# DESeq2 is loaded\n",
    "\n",
    "library(\"DESeq2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# Builds a sampleTable based on the data files in the tumor_samples directory\n",
    "\n",
    "rna_data_directory <- 'data/rna/tumor_samples'\n",
    "\n",
    "sampleFiles <- list.files(rna_data_directory)\n",
    "sampleNames <- gsub(\".tsv\", \"\", sampleFiles)\n",
    "sampleCondition <- c(rep(\"tumor\", length(sampleNames)))\n",
    "sampleTable <- data.frame(sampleName = sampleNames,\n",
    "                          fileName = sampleFiles,\n",
    "                          condition = sampleCondition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# A DESeq dataset is constructed using the sampleTable and the path to the tumor_samples directory\n",
    "# The design is set to 1 to indicate that there are no factors dividing these groups (this could be changed if desired)\n",
    "\n",
    "dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,\n",
    "                                       directory = rna_data_directory,\n",
    "                                       design= ~ 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# The DESeq dataset is transformed using the log-like variance stabilizing transformation (http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)\n",
    "\n",
    "transformed_dds <- vst(dds, blind=TRUE)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# The transformed_dds_matrix object is a summarized experiement object\n",
    "# This is just a matrix of the transformed counts\n",
    "# This matrix is then saved to the temp folder for passing off to Python\n",
    "\n",
    "transformed_dds_matrix <- assay(transformed_dds)\n",
    "write.csv(transformed_dds_matrix, \"data/temp/transformed_dds_matrix.csv\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The transformed_dds is loaded into the Python environment with Pandas\n",
    "\n",
    "transformed_dds_matrix = pd.read_csv(\"data/temp/transformed_dds_matrix.csv\")\n",
    "transformed_dds_matrix = transformed_dds_matrix.rename(index=str, columns={\"Unnamed: 0\": \"Ensembl_ID\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The \"high expression\" and \"low expression\" cutoffs are calculated for each gene\n",
    "# Here these cutoffs are based on the median expression level for each gene \n",
    "# But can be based on other thresholds (i.e. upper tertile/lower tertile)\n",
    "\n",
    "sample_columns = transformed_dds_matrix.columns.values[1:]\n",
    "high_expression_cutoff = transformed_dds_matrix[sample_columns].apply(lambda x: np.percentile(x, 50), axis=1)\n",
    "low_expression_cutoff = transformed_dds_matrix[sample_columns].apply(lambda x: np.percentile(x, 50), axis=1)\n",
    "\n",
    "transformed_dds_matrix['high_expression_cutoff'] = high_expression_cutoff\n",
    "transformed_dds_matrix['low_expression_cutoff'] = low_expression_cutoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download clinical data from http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/PAAD/20160128/gdac.broadinstitute.org_PAAD.Merge_Clinical.Level_1.2016012800.0.0.tar.gz\n",
    "# File should be extracted and \"GBM.merged_only_clinical_clin_format.txt\" should be placed in the \"data/clinical\" folder\n",
    "# The clinical file is read in with pandas \n",
    "# The dataframe is transposed to force the clinical feature titles as column names\n",
    "\n",
    "clinical_dir = 'data/clinical/PAAD.merged_only_clinical_clin_format.txt'\n",
    "clinical_df = pd.read_csv(clinical_dir, sep='\\t', index_col='V1').T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Patient barcodes need to be converted to uppercase for later merging\n",
    "\n",
    "clinical_df['patient.bcr_patient_barcode'] = clinical_df['patient.bcr_patient_barcode'].str.upper()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The most recent recording of 'days_to_death' is assigned to the 'death_days' column\n",
    "# This should be an unambiguous number, but the max is taken anyway just to pick one value if there are multiple entries\n",
    "\n",
    "death_df = pd.DataFrame()\n",
    "death_df['patient.bcr_patient_barcode']  = clinical_df['patient.bcr_patient_barcode']\n",
    "for column in clinical_df:\n",
    "    if 'days_to_death' in column:\n",
    "        death_df[column] = clinical_df[column]\n",
    "        death_df[column] = death_df[column].astype(float)\n",
    "        death_df['death_days'] = death_df.max(axis=1)\n",
    "\n",
    "death_df = death_df[['patient.bcr_patient_barcode', 'death_days']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The most recent number for \"days to last followup\" is assigned to the 'followup_days' column\n",
    "# Days_to_last_followup is defined as the time interval from the date of last followup to the date of initial pathologic diagnosis\n",
    "# This means the most recent followup would be the largest number, which is why the max number is taken\n",
    "\n",
    "followup_df = pd.DataFrame()\n",
    "followup_df['patient.bcr_patient_barcode']  = clinical_df['patient.bcr_patient_barcode']\n",
    "for column in clinical_df:\n",
    "    if 'days_to_last_followup' in column:\n",
    "        followup_df[column] = clinical_df[column]\n",
    "        followup_df[column] = followup_df[column].astype(float)\n",
    "        followup_df['followup_days'] = followup_df.max(axis=1)\n",
    "\n",
    "followup_df = followup_df[['patient.bcr_patient_barcode', 'followup_days']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The recently created death_df and followup_df dataframes are merged\n",
    "\n",
    "df_list = [death_df, followup_df]\n",
    "all_clinical_df = pd.DataFrame()\n",
    "all_clinical_df['patient.bcr_patient_barcode'] = clinical_df['patient.bcr_patient_barcode']\n",
    "for df in df_list:\n",
    "    all_clinical_df = all_clinical_df.merge(df, on='patient.bcr_patient_barcode')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the patient did die, the 'death_days' value is recorded in a column labeled 'new_death'\n",
    "# If the patient did not die, the number of days since last followup ('followup_days') will be recorded in this column\n",
    "\n",
    "all_clinical_df['new_death']= np.where(all_clinical_df['death_days'].isnull(),all_clinical_df['followup_days'] , all_clinical_df['death_days'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/brett/.local/lib/python3.6/site-packages/ipykernel_launcher.py:5: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \"\"\"\n"
     ]
    }
   ],
   "source": [
    "# An alive_dead_df dataframe is created to record patient that lived/died\n",
    "# This dataframe has a column named 'death_event' which is coded '1' for patients that died and '0' for patients that remained living\n",
    "\n",
    "alive_dead_df = clinical_df[['patient.bcr_patient_barcode', 'patient.vital_status']]\n",
    "alive_dead_df['death_event'] = np.where(alive_dead_df['patient.vital_status'] == 'dead', 1, 0)\n",
    "alive_dead_df = alive_dead_df[['patient.bcr_patient_barcode', 'death_event']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The alive_dead_df is merged with all_clinical_df on the 'patient.bcr_patient_barcode' column\n",
    "\n",
    "all_clinical_df = all_clinical_df.merge(alive_dead_df, on='patient.bcr_patient_barcode')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The transformed_dds_matrix is subset to only include a single row corresponding to some gene of interest\n",
    "# In this case the gene is the c-Met receptor \n",
    "\n",
    "gene_id = 'ENSG00000164078.11'\n",
    "gene_subset = transformed_dds_matrix[transformed_dds_matrix['Ensembl_ID'] == gene_id]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The gene_subset dataframe is transposed to make the samples rows and a new header is created\n",
    "\n",
    "gene_subset_transposed = gene_subset.T.reset_index()\n",
    "new_header = ['patient.bcr_patient_barcode', 'normalized_counts']\n",
    "gene_subset_transposed = gene_subset_transposed[1:]\n",
    "gene_subset_transposed.columns = new_header\n",
    "\n",
    "# Patient barcodes are truncated to ensure correct merging with clinical data\n",
    "\n",
    "gene_subset_transposed['patient.bcr_patient_barcode'].iloc[:-2] = gene_subset_transposed['patient.bcr_patient_barcode'].str[:-4]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The high and low expression cutoffs are obtained from the gene_subset_transposed dataframe\n",
    "\n",
    "high_expression_cutoff = gene_subset_transposed.iloc[-2, 1]\n",
    "low_expression_cutoff = gene_subset_transposed.iloc[-1, 1]\n",
    "\n",
    "# These entries are then dropped from the dataframe, because they are no longer needed\n",
    "\n",
    "gene_subset_transposed = gene_subset_transposed.iloc[:-2]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gene_subset_transposed is merged with all_clinical df on the 'patient.bcr_patient_barcode' column\n",
    "\n",
    "all_clinical_df = all_clinical_df.merge(gene_subset_transposed, on='patient.bcr_patient_barcode')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A new column 'event_rna' is created in the all_clinical_df dataframe\n",
    "all_clinical_df['event_rna'] = pd.Series()\n",
    "\n",
    "# Every sample is coded as expressing the gene of interest above or below the median expression level\n",
    "\n",
    "all_clinical_df['event_rna'] = np.where(all_clinical_df['normalized_counts'] > high_expression_cutoff,'high_{}'.format(gene_id),all_clinical_df['event_rna'])\n",
    "all_clinical_df['event_rna'] = np.where(all_clinical_df['normalized_counts'] < low_expression_cutoff, 'low_{}'.format(gene_id), all_clinical_df['event_rna'])\n",
    "all_clinical_df['event_rna'] = all_clinical_df['event_rna'].replace('nan', 'medium_{}'.format(gene_id))\n",
    "all_clinical_df = all_clinical_df[all_clinical_df['event_rna'] != 'medium_{}'.format(gene_id)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Any files that are currently in the data/temp folder should be deleted\n",
    "\n",
    "tumor_sample_dir = 'data/temp/'\n",
    "files = glob.glob(tumor_sample_dir + '*')\n",
    "for f in files:\n",
    "    os.remove(f)\n",
    "    \n",
    "# The all_clincal_df is saved to the temp folder\n",
    "# The name of the gene that is being analyzed in incorperated into the file name as well\n",
    "\n",
    "all_clinical_df.to_csv(\"data/temp/{}_all_clinical_df.csv\".format(gene_id), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: ggplot2\n",
      "\n",
      "R[write to console]: Loading required package: ggpubr\n",
      "\n",
      "R[write to console]: Loading required package: magrittr\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# Load the survival library for analysis and survminer library for plotting\n",
    "\n",
    "library(survival)\n",
    "library(survminer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# The name of the all_clinical_df csv saved in the temp folder is gotten\n",
    "# This file name is subset ot get the gene name\n",
    "\n",
    "file <- list.files('data/temp/', pattern=\"\\\\.csv$\")\n",
    "gene <- strsplit(file, '_')[[1]][1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] \"data/temp/ENSG00000164078.11_all_clinical_df.csv\"\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# The full path to the all_clinical_df csv is constructed\n",
    "\n",
    "all_clinical_df_path <- paste('data/temp/', gene, \"_all_clinical_df.csv\", sep=\"\")\n",
    "all_clinical_df_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# The all_clinical_df csv is read in\n",
    "\n",
    "all_clinical_df <- read.csv(file=all_clinical_df_path, header=TRUE, sep=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Call: survfit(formula = survival_object ~ event_rna, data = all_clinical_df)\n",
       "\n",
       "                                   n events median 0.95LCL 0.95UCL\n",
       "event_rna=high_ENSG00000164078.11 89     34   1059     634      NA\n",
       "event_rna=low_ENSG00000164078.11  89     24   1502     738      NA\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# A survival object is created by indicating the time to death/last followup and death event data \n",
    "survival_object <- Surv(time=all_clinical_df$new_death, event=all_clinical_df$death_event)\n",
    "\n",
    "# The survival_object is used to fit the Kaplan-Meier curve with the survfit function\n",
    "\n",
    "fit <- survfit(formula = survival_object ~ event_rna, data = all_clinical_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nOzdeVxU5f448M85ZxaGfVPcAUHFfYlyzRUVU69LZYalluutNDNLU283+1bmVX9l2KJZ6XVpVTMtNbmWRqioiBsuKC6ICozAwAzMzDnneX5/HJ1oGGCAw8wBPu/Xfd1mnjnnM5/nYfhwPHPO8zCUUkAIIaQ8rLsTQAgh5BgWaIQQUigs0AghpFBYoBFCSKFU7k6goduyZUt+fr6TG4uiyDAMy8r2Z5UQIlc0SqkgCGq1WpZoACCKIsdxMkaTcegopZRSGaOJoqhSyfbLKO/QCYLAsqyMnQUAhmHkiqbwoeM4zvnORkdH9+7du3QLFmg3y8/PnzNnjpMbFxUVqVQqnU4ny1uLosjzvIeHh1zRDAZDYGCgLNEAwGQyeXl5yRWtsLBQo9HI2FlBELRarSzRBEEoKioKCAiQJRrUwtBptVoZO0sI0Wg0skTjed5oNCp26AwGg06nc7KzhJDPPvsMCzRC9ZyMJQYAfH19ZYwm49EuAKjVahmrM8g9dDWH56ARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaIYQUCgs0QggpFBZohOqVn376qVevXn369OnSpcvkyZONRiMAbNu2zZl9ndwMuQwWaITqj5KSkueee27nzp1JSUmnT5/29vbesmVLcXHx6tWrS29GCCm7b9nNkNthga5LMu/d234p/fcbmdfNlowSs7vTQYpjsVgsFosoigDAMMwnn3wye/bsWbNmpaenT5ky5eDBgyNHjhw8ePCGDRtEUZwxY0bfvn179+49d+5cALBtVvYl5C54J6EMCCE//vjjzp07N2/eDADFxcXx8fEMw2g0mjlz5lgsltJPa3Kn/7fp6ZuBbZmfPxRAoPStsFbydQLVB/7+/m+//XbXrl379OkzZMiQcePGhYaGLlmy5OzZs5s2bUpMTDxy5MiVK1cCAwOzs7OjoqI+//xzAOjWrVtqaqpts7IvdevWzd09a6DwCFoGBoOhXbt2zZo1k54eOHCgZ8+er7/+eosWLZKTk+2eujdVVO/Nnz8/IyNj2rRp6enpPXr02LJlS+lXO3ToIM2XEhQUdO/evXHjxk2cODE7O1uv19u2qeAl5GJ4BC2DgICA0hMC3Lx5c+zYsQAQGRl54cIFvV5f+qk0Gcqvv/566tQpABAEoaSkxMk3svB8Cae5Ybb869pNAFh2PbP0qwzAk0GBGyPDnIxGCCGEyLXmmRTK+b5UqkojUylpZigZOyuRK1r1hq7stFmEkMLCQn9//7Fjx44dO3bEiBHLly//4osvbBvYpovauHFjamrq7t27OY7r06dP6SAVvIScV1JSQgixWq3SGadKOfxEYYGuRYSQ0jMNln6q1+szMjIAoFmzZpX+/PaeP/+s0WrhGPD0BfgrIAuUlHpKAX4tMDj5aQAASikhxPntK40GAHJFgwczScoYTd7OyhhNKtCyRNu7d++iRYsOHTokHSanpqZGRkayLGuxWOy2zMnJiYiI4Dju1KlTFy5csFqtts3KvlTzxBogURSlz4mT2zs8WsICLb/w8PD09PSWLVump6eHh4d7e3uXfiptExcXFxcXBwDx8fHe3t4VBxzZtWtSTu4Vo/GbG7cOqVT+PN/a24vT6V5s0dxPpwtUcZ4cd85UPObcBYZhKo1mI/t0ozzPO//ulWpo043KMnQjR47MyMgYOHAgx3GEkK5du3744Yeenp4MwwwaNOj//u//bFvGxcWNHz8+JiamU6dOixcvXrRo0S+//CJt9tVXX9m99MgjjwQHB9c8vQbF29vbYDB4eHg4P91o2S+oGFzVu+YuXry4ffv2tLS0Dh06DB48uFu3bh999BEhxMvL66WXXrJYLKWf2k18Hh8f7/x80P8+cGCzQFqWFA/mGOLju2zQIHhwSH6yyBh98rQHyy5s1XxO82ZB6sr/9OJ80NWm8PmgkUJUYz7oF154oXQjHkHLICoqasmSJaVbFi5caHus0+lKP60xBgBozl2Sm03btWGat5RaNSwLAGZCll3P1DDs4tAW8r0jQsg98CqOumRBr17nw1v8z1P1xo0rS9PPU/Gv01udvDxXtA5r66kDAAuV58srhJB7YYGua0LDVbGjGekfTYa/FjNkAF5v1byPrw8AFAqiiGeuEKr7sEDXTQwDAPzWr2j2nbIvfnjr9uPnL7o8J4SQzLBA10nsQz2B44BSWlhYur2Pn6/0Ez1jLHZLYgghGWGBrpNUj/2D8fMv2z6jacivXTu6Ph+EUG3AqzjqNv6/Gxi1in2op2rkWHfnghCSGR5B11WMfyAAgNVCTSaSepLm3aP5eTTvnrvzQgjJBo+g6yr11JlUn0OyMoXt34DVKp5MZliWEgLRvd2dGkJIHlig6yytlmnekrFYAIDyPLlyiQGGAtDwSACwUHKyyNhYo24p091uCCHXw1McdZx0874o0OsZ5PpVev0qOfInANy2WKNPnm5zLOWu1QoAOMF/Q4bnvmSRmJhod8MwALz00ktnzpypvTfFAl2HkbSz/JefAccBMA/+B/3Pp+oTdk7Juc0AWAjJsfIAsCU757/ZOe7OF7kHSTkunnR2IvJ9+/ZFRESMeuDo0aP79u0LCQm5d+8eANy6deuZZ54xm82TJ0+eMGHC448/PnnyZJ7nAWDTpk0xMTGxsbH9+vXbs2cPABgMhieffPKpp56aOnWqIAg1eSqK4n/+85/GjRtLSVJK58yZ88QTT4wcObLwwZWmW7Zskebytdt348aN0uSr3bt3/+STT65cuTJy5Mjnn39+7Nixubm50r6l49u9l82yZctsBfrbb7+VHrz99ttlq7aMsEDXYWybdpp5CzWvLOIe6c093JOJiAQAtnFI4JwFG8eO8arB0i2oIRs3btyeB3r16gUAgwYNeu2112wbJCUlBQYGfvfdd9u3bx81apRerz927Ni2bdt++eWXffv27du3T1qY4ssvvxw7duy3334bFRW1e/fumjzNycnp3bt327ZtpQQSEhL8/f1/+OGH999/X/rzkJ2dnZCQIM2uZbfv1KlTf/zxx+3bt0dGRk6dOvXgwYPjx4//8ssvu3Tpcvz4cSlg6fh27yU5e/ZseHi4p6en9HTv3r3Sg8DAQJ1Ol5n5t5nZZYTnoOsytYYJCAQAxj+AEsKYzRSAZt+xfrxa8+r9v+pPX7isY9k7VmugSrUktKW61PzUqF4SvttCrl39W1NJCQUgKfYH0eoX5zPePmUj7Nq168qVK9LjrVu3AkB0dPSNGzcSEhKioqIA4OGHH169evXLL788YMCA4cOH+/n5bdiwYerUqdK0bd7e3m+//TYAnD9/PjY2FgAeeuihpKSkzMzMaj8dN25c06ZNbRkmJyebTKZZs2apVCppods33njj3XffnTJlStn3HTduHACsW7cuLi7O09MzLi5uyJAhe/fuNRqNb775phSwadOmtvilH9ukpKT06NEDAL766qv9+/cfO3Zs4sSJPXr0eP3113v06JGSktKyZcuq/JSchQW6/mB8fIFTgSiAxULz9EFqlVEU00z3bym8bbEGJR5lGWZ928gJjXFu33qLFhU5PONMS8rcXFrOXPJjxoxZtWqVXeN77703atQoaaFCHx+fn3/++ebNm4cPHx45cuTatWsFQZBmpj9w4MDq1av1er20PueD9yl35YqqPpVYrdawsLC5c+d+8MEHmzdv9vPz69WrV/Pmzcv07/6+hJCvv/768OHDALBmzZp58+Y9/fTTq1ev3rp1q1TTK1VQUBASEgIAzz333HPPPTd16tSNGzdKL/n5+RUUFDgTpBqwQNcH3JDh0gPVmCes771JDQUAcLBrp6TCokJRBIC56RkipT90bN/WU9fIiamiUd2lGjGa9h8MACAI1FwCAOT8GaCU7dQVABgPHajufwAYzypMt+3j4zN//vx3330XAPbt26fT6QYMGPDMM8/k5uYePXp05MiRCxYsGDdu3NChQ4cOHfrss89aLJauXbseP368ffv2ycnJXbt2DQgIqPZTu2Q6deoknVXw9fU1m80///wzwzBJSUlpaWnr16+3e18AOH78ePfu3aVinZeXFx0dDQABAQF5eXlOdt/Pz89gMDh8yWAw2NYjlR3+rtZTVr61zqO17v789PPSM0SAc6ZigyhoGXZooL+Oxa8f6iemWQv7A868e5QQrnu0kxF27tx58eL9ybaeffZZPz8/6fGYMWM2b97s4eHRvXv3F154Yc2aNQzDaLXaTz/91M/Pb+bMmbGxsb6+vhzHDR48uEuXLhEREdOmTdu1a5efn9/ixYuLi4ur/fTIkSMrVqy4cOHC2LFjJ0+ePHbs2BkzZkycOLG4uHjz5s0vvviilGFMTMzMmTOLiopK7wsA586di4iIkLZ56aWX5s2b9/333xcUFKxfv/7EiRMffvjhiy++aIvfpk2b9PR023uNHz8eALp167Zu3TrbENkOnwEgJSVFWh2pNuCKKm5WpRVVioqKVCpV2aVCS7Oueofm5nADYlSP/cPWqDmcxJO/ftCvtmy+KiIMV1Sptrq1oor4v/2UENXQEbUUv4GIiYnZs2eP3UcoLy9vypQpu3fvdrhLzVdUwcOo+oaNbAcAQP62Aun8Fs0f8vF+yMe7qUYDALk875bckFtwQ4Zjda65N99885133nGmUUZ4iqPekb5R+ftKzO+3DgUIBYC3rt9cdr22LglCqB7r379///797RrXrl1bq2+KR9D1k3jsT3L5QnmvSnevIIQUDgt0fcNGdQCNBgBoQX7ZV0M0GgBI+vs0/wghZcICXd+w7TqwrcLKe1VatBC/GEaoTsACXW+JB38VE/Y6fIkACFikEVI8LND1kacXAND8POH3/9m9wjIMAJhE8Znzl3F+O4QUDgt0PaQaN0E1ejyAg3MZ7XS6zl6eALAr/96WXL3rc0MIOQ8LdD3EeHqxbaIcvqRhmQ5eni7OByFUPVig6ymWAQAQBZqT7fB1C6Hxt+92SD0XcexkxLGTz1647NL0kAvhcg2ysJuwv7an6pfgjSr1E+PtCwBAqXjimO2e7wnnLn2vv39ag1LI54V8XpCeZpSYY/z9pjQNcUeyqHZtyc4RKH0rrJUzG+/bt+/FF19s37699HTp0qUFBQVTpkxJS0sLCgq6devWokWLNmzYMHPmTLPZLIqil5fXF198oVarN23atHnzZpVKZTQaFy1aNGrUKIPBMH36dJZldTrdhg0bTCZTtZ8yDLN69epVq1bl5OQAAKV07ty5d+7cKSkp+frrr319fQFgy5YtH3/88ZEjR+zed8uWLT/++CMA3LhxY8aMGcOGDXv55ZdDQkLy8vI+//zzRo0aAYAoirb4pR+XHplly5bt2rXL9vTtt9+u4CZvueARdD3laNrn1RHhW9pHfRQZ4a9SaVn2ER+fF5qGLAu7P4/tmzfwDkMEgBP2OzFhP9T+VP0SPIKu58jF8/DgCLqlTjNJFwQAH2bdyuNpbKDfG01D8ln239czAYCnhDou7KguWXPrdlpxSemWlCIjBbhT5vbRFa1D/VUOKgBO2F/BhP2l1epU/RIs0PUQSTvLf7sFGAYopdl3LG++Dgyjih3FtuvABAYBAAhC2cna71j45TduLQ5t4YaMkXz25uXvz3Mwf/zJIqNdy79DW/o7KgA4YX9Ztgn7P/zww3379j311FPPPfdcrU7VL8ECXQ9JaxXSnGz+y08BgHuoJ9d/EDl+VDyZLM1qdjHtJCXkl6i2vkdTePjrUrwl124suXajp7fPqjbhIiGhOo/kQmOIWsWxbGqRqZFa9bCv920rn2W2MAwTolaF6jwIpbZZp5ESTGgU3M3bGwB4QosJAYAUo5FQiPbxBgBPllE/mArcuyqrVuKE/dKE/fPmzZs3b57UWKtT9UtwPmg3++CDD2bOnOnkxlarlWVZlaN/ljpgKFCtWQEAtE0UadueTb9EKaVtowCAO/onpbRw9PjTfsHnrZZfDUUJhkKeEl+Oa65Re3OsF8uKAAN9vDfr81pqNRxAanGJH8dNCg5MKjJes1gZgJZazUAfbxFgcbMmAEAptVgsck24DAA8z6vVarmiWSwWjuOcHbrKUEoJIZxMy/ISQnier8bs0k7Ol/3Ojcxqf0koTdh/7ty5BQsWAMATTzzh4eGxevXqF154gVJaesL+LVu2rF+/3jZh/5w5c0wm07Rp00RR9PPz+/zzz6U5+Kv3NDk5ecWKFX/88cejjz46efLk0aNHz5gxw2w2SxP225YUiImJSUhIkCbst+3LcdwXX3xhNBpffvllALh27dq8efMaNWokTdifkZFhm7Bfii9N2G97L2nC/tTU1HXr1n366aelx2rChAmrVq1q1crxwJpMJovFolKpnPycEEI2btxoNzs8Fmg3k33CfhtaaLC++6+KtlCp6eK3Pby83rmR+e6NTDOhgWpVqFYrAuWAEYH28Pben5ffSKPmgEkvKfHmuKGB/ueMxXesVgagkUbdw9vbg2U/bRsBOGF/DdT2hP1VKtCoPHYT9lc8Vb+k5hP24ymO+un+aWiVCkQR6IMv/xgKpf8cCzzz9mLx0YHQ+f5iSHm8kPfgwjsAOGMshlJfLhlFsvluru3VO1Ze2uDlFs2iPJ36m4HcYmloLX6L1XBIc/Pbpuev7an6JVig6yfpNDQAiIm/gcXCRvciJ5OBiky7jgAg/m8fvXMbAOiU6f/mdO9du1mT9+py/JQnx34W2XqYBj9OqN6ym7C/tqfql+BvVD2l1jABgQDAeHpTD082rDW9mk4J4Tp3AwCSeEg6kmY8dPNahMU2CTEI5NWr1wRCBSACoSwwPNBwD4+zRpMnx7EM5Fl5Nce21HjcE3gzEVkGfDhOz/NmQj5pGxET4B/IMoLR/joBhFBNYIFukFgGgAGg9PiR4NaR/TQ+APDBLZVAqYphbP8fE+B3lze30mpVDHOSCP4q7qmQgERD4dUSM8NAK602TxSAQFudLsxDK4qi41XpEULVhQW6nuOGDLd7AADqWXOFHd+Ix5KYcmbqcJJVFCvfCCFUXVigGyq1/TfLCV07ld2q4u+XPA4nCZSeLy7u7+8rZ24IIQDAuTgaupxsKC6u9t7S9P+LMq7Llg9CqBQ8gm5w7l+BJ52dsFos770JLKt+4mm2S3cnIyzNuPnezUyA+9fsFQoi+/ufc5s1eTO4tq7kRahhwgLd4NiuwLP+5/+AiFynLuDlzUZ1dD7CvBbNYoMCAOCfl6+eM5lULPNb186RGhWUlFS6L0LIeVigG54HV+BJqNkMDEtuXAOGYVuFgRN3PQVrVNKFHz4qFgAYCv38fERRNGCBRkhWWKAbLgaAApAL5wCApCQDABfdU/XkpKrGIYCzBSBUK/BLwgaHpJ21/HuhZekC+rcZRxkAoLezqhTqy3ZtAIAAfHTrTrFoP38pcheTySRjtMLCQovFIlc0QRCsVqtc0Xiez8/PlysayD10NYdH0A3O/XPQAi8cPkgEkVOpgBLw9Rf/tw+8vKsUSs0wAEApvHwlQwV0oifOO4qQnLBANzy2c9C+/lQQGI2GEsJ4ewMALTKQs6ls526VxvhJnzf54mUAUDGMSIECfTXjukeLZs/LN5sdQghPcSAAAGBZAKB37/BbviTXMyrdPCbAPzW6e2p095dbNJMm7J/fvPkw36odgCOEKoZH0A0XO3iYyPOchwcA0KJC9splcuE88FaamwthrSve15NjwzgtAPirODUDABCg4jxZ/HuPkJzwNwoBADA+vupJzzFqDQCQK5fcnQ5CCAALNCqL3skS9v5EzqY6s/HS0JZRnp4AIOLFdgjJDU9xoFJYBgBo9h0x+47IcdqOXcCJsxbSksuf3Mme4VPJiRGEUJVggUZA0s7y2zYBz/+tVRTFI3+w7TsBwwClTGBQebuPCQ7cqb+XbeXL2wAhVD1YoBGwbdppXn6NmkzkZDLJymSCG5HTKQAARUXiyWSGZSkhqqEjytv9EV8f1+WKUEOCBRoBqDVMoxCmEdCr6cy9XMbn/iL25NoVCsAAQ4GSFi3Z9g4mjLYhFK5brIVmC6FUuvAOIVRDWKBRGQ/OO0sXREvf/Qnfb9W8ubyCnUQg3xYYfKy8QOlbYa1qPUmEGgC8igP9hRsyXD1zDhveGtQaYFkA5sH/gJqKLf96jZw5VcHuZkpLCDETki8IeEkHQjWHR9DIHtumnebVxSDw4uHfgBJaUkLOnQaGYZo0rWDaaELho5x70uMVN7OebxryRbtIV6WMUP2EBRqV8WCyDsY/gBLCqNT3280WUKtLb/iTPu/ptEslxME8dpvuZg8P8J/QOLj200Wo3sICjSrDMAAAlNCcO+RqOhvZ1vZKTID/iYe63ROEL+5kHy0sClerdBrNXas1yVDYwdNrZBCugIVQjWCBRuXihgwX9+8Rkv54MLk/8J9/DADcw71UTzwNAJ4c295LBwC/FxTcMJd01mp8PD3PmCgAaFjGi+Pcmj5CdR4WaFQR7tGBTLsOAMCviwcicrGjIO8eeOjcnRdCDQIWaFQhT282zBsAGIahAIxOR7VaEEWwWkuvXrg0tOUbLZoZDIbAwMDvc/XfP/i2ECFUE1igURUIO7+THpAzKZpFb4Gq3M9PviB8n6sHABXDxAYG6HAmUoSqDgs0qghJO8t/uwUAgIil22lRoeXtxeonnma7dLfbhQUGADJKzBPO35+29M2wlsvw1hWEqg4LtPz+97//HT16FAByc3OHDRsWHBy8efPmJk2aAMDMmTMbNWrk7gSr4P4ChgBi4m9gsbDRvciJo+LxowCgnvkS27hJ2V0G+PuODAq4a+UB4JbFkm3l83jBxWkjVD9ggZbfkCFDhgwZQghZuXLlkCFDkpKSxowZExMT4+68qsV2TbSnN/XwZMNa08sXpVcYP//Sp6FtgtXqPZ07SI9fvXrt/2XedlmyCNUzWKBry759+/r376/Vak0m04kTJ06dOuXj4zNt2jS1Wg0AiYmJaWlpAGA0Gs1ms5MxRVGklDLShck1RggRRbHy7QAAgBEESohgNjPC/cNh6/fbqEoFLUOh96NSNEqpXV8EQQCAAiu/NetuIxXX368K894JguD8yFRKFEWel21CVEKI1F+5opUdupqQd+gIITzPy9hZiSzRpN8IhQ+dk511OCxYoGsFpfTw4cPLly8HgN69e/fq1Ss4OHjHjh0JCQkjRowAgIsXLyYkJABAhw4dnC8c0s9PrkJDKa3C70mf/gAAPA99+msSfwNRhEtpDACknbV2fxg4TvoFtstNin+x2JRYWNhCo/7FM8L59KQPt/PbV4xS6vxfI2eiSWSMJu/fD3mjSXVQlmjyDp30t02xQ0cpFQTB+QJddliwQNeK9PT01q1bS4e6WVlZQUFBAKDT6Ww/qunTp0+fPh0A4uPjfXycPa4sKipSqVQ6nTyXIUsHlR4eVZ4alDw3m2ZlUotZPPgrUOqlz2Yj24miaDAY7Pqi0WgAgONUDCEcxznfUwAwmUxeXl5Vza08hYWFGo2mGp11SBRFQRC0Wq0s0QRBKCoqqtLgVEz2odNqtTJ2lhCicXRyrBp4njcajYodOoPBoNPpnOwsIYQrc28XFuhacePGjaZNm0qPfX19165d6+PjwzDMvHnz3JuYLNg27aBNO1poEA/+CgBCwn5NZLsKtr9mNpcQYiE0udD4iK+3q9JEqM7DAl0rhg4danscFhYmneuot6wWcPQv1p/0eWsy7wBAjpUHgCJB7JlymgF4uXmzD9qEuzpJhOogvH0AVdeD7yppVqb424Gyr8cE+P/atdOExo3ae3lpWBYAunh77+/S8Y3Q5i7NE6E6Cws0qg5x/x7rO/96sKI3CPt/FhbP1x7cX3obT44dHOD7bYe2cY2DfDkWAJpoVEMD/RvLdP4RoXoPT3Gg6vhrEqWtX0GhgWvfEfoPLtHgUoQIyQkLNKoW2yRKHEcBwMeHCQ0nBoPDbZeGtgxWq/95+WoJIRklZlxSFiEn4SkOJAfnLjG+Y7H+NzuntnNBqN7AAo1qhH2oJwCQK5ed2dhCyB2r9YzRJMh0nwJC9RsWaFQj91fAspjptasVbCZ9mZhpsa6/nd31ROr0S1dckRxCdRwWaFQjjEoFANRiJQn7ytvmJ33egqvXWAZsc4hsyc75LkfvkgQRqsPwS0JU62IC/M8+3AMA1mbdOWcy7c8r6O7tg0vKIlQpLNBIDpTS/LzyXvTk2DBOCwD+Ks6TYwFAi0vKIuQELNComkjaWX7LVyBKs49SKMjzWvGWBRjVowO5UePcnBxC9QIWaFRNbJt2mjnziV5PEn8n1zOAU1mGDPcODWebNCtvl6WhLTdn5+zMzbMQihdEI1Qp/JIQVZdawzRtwXXuxgQ1AgCglDGZ2Mi24F35fHU5PF4QjVDlsEAjNygUhD8Nhetv3+XxgmiEyuegQN+4cWP+/PmdO3f29/f39/fv1KnT/Pnzb9y44frkUJ3Ado8GAKBUffKYdeU75FJaBRurGAYACgQxId8w6/LVH3LxYjuEymVfoNevX//YY4+Fh4dv3749Ozs7Ozv7xx9/bNOmzahRo9avX++WFJGSkbSz/OYvAAAoYcwlVJ/Db9pAzpxyuPFP+rzZl69qWIZjGOmq6GmXruAF0QiVx/5LwkuXLqWmpkoLm0oiIyMjIyNnzJixcOFC1+aG6gC2TTvNK4uoPlc4+qd4+xabp+eaNQcPHc3PA0qZwCAAoPl5tCCf8fOPCQg4Hd0dANZm3fkuV59ptrwZ2goviEaoPPYFevXq1dIDSmlSUtK1a9dCQkIGDx6sUqlsLyH0F7WGCQhkAgKZ6xm0IA/yQLydRX47wEW2pYSoho4AAJJynFy9zIRHeg4dYbsgWsUAAASpVXhBNELlKfdLwqVLl27fvr2goOCPP/6YPHmyK3NCCCEEZY+g33nnnddff12j0dy+ffvzzz9XqVQA8Mwzz0YtwPsAACAASURBVLgjN1SXsIOHEWMRd+smiCK9kyUa8oECfyMDAGjePSgupvn5wj296vGnQK1ZGtry9wLDtRKLu7NGSNHsC3R4ePjw4cOXLVs2YcKEgQMHUkoZhomLi3NLcqiuEPfvEQ7++uCLCwolxbSkGABoXqkvAEuKxTy9eDZV/dQzbJfu7kgToTrGvkBPmjTpsccee+ONNwRB2LVrV1BQkFvSQnUL9+hA2jK0JCtLm5xITcWg0zGeXsBxbGg4UKB3b1NDATUWgcCrxj/FRnW07fhHQSF9cCl0L1+fLt5ebuoBQkrk4FbvgICAzz777MiRIxMmTHj22WenTp3q8qxQXePpzbTrwDdp7nH9KtzKhMZNSn9JKP5vP7l6mZpLQOAZ/wDQaACAYxgA2Jyds/nBLYWROo/0ng+5sRMIKY39l4SHDh0aN27co48+un79+nXr1un1+pEjR164cMEtyaF65e83Db7aovmIwICYAP+YAP+HfbwBoMi5dbMQajjsj6A3bNiwYcOGwMDAS5cuLVu2bPPmzRMmTFiwYMF3333nlvxQ3aKe8WLZRm7IcG7IcOu/F5au0MMC/YcF+kuPjxQW9Uk545IEEapL7As0z/OnTp3y9/dPT0/38PAAgFatWmF1RrIRBHdngFCdYV+g33vvvU2bNun1+oiICLwzBcmJYQCAnD7Jtmvv7lQQqhvsC/Tdu3eXLVvmcNOkpKQ+ffrUfkqoXiFpZ/lvtwAAWC0AIKaeFM+fVT8+seyVdiKFwwWFLbQanCcaIYn9l4TLly9fuHBhdnZ26cbc3NxFixa9//77LkwM1RNsm3aaeQs18xayYeG2p6WvtLMxE/GtGzdxnmiEbOyPoHft2vXJJ5/07NnT39+/VatWAJCZmWkwGF599dX33nvPHRmiOk6tYQICAQBU6r89LYOn9LbFcrLIuP72XR8V95inzoVZIqRE9gWaZdmXXnrppZdeunjx4s2bNxmGadWqVbt27dySHGogOGAAwELopWLzpWLznnv5APB561bTnVicBaF6rNw1CaOioqKiolyZCqrf2Mi25PLFsu0/6fOeTrtUtv1woXF67WeFkJLhorFuJgiCyWRyfmNCCCFElreWQoky3R5CKaWUVtAX1mplAUjGFfO6j0CtIUNH0MBgAOitUf/Rod09Qdh2Lz+1uLizTnfJbEk2Gttp1c6PTKUEQaCUythZQogg0yWDhJCKh66qeJ6XMZooihaLRcbOAgDP83JFI4QoeejMZrOTnXX4y4gF2s1UKpWXl7MTUBBCVCqVTifPyVlRFHmel652lyWa1WqtoC+ir58AACYjk3EFALStQrmYEQDgBdAIAACOWSw3BL61t+c9SsAIOpXa+ZFxJj2NRiNjZwVB0Gq1skQTBEEQBBk7azKZ5B06rVYrY2cJIRqNRpZoPM8reegEQfDw8HCys4QQrszc6I4L9C+//BIbG8uyuKQskg33UE/QaMDKk5PHyI1rFJeLRagyjktwfHx8ZGTkv/71r2vXrrk4IVRvqVRc94e5nn2YRo0dvr40tOXBrp3eCmtV7Xe4brZklJhLP7a13LTy18yWspuV3vdwQWHZ9vLcsFgdBnE+AkKVclyg9+7de+LEidDQ0NmzZ8fGxn7zzTdynTNCyEY8myok7BdPn6J592QJuCU7x3YZtfTY1vJtXv5W/b2ym5Xet0pXYW/L1TsMgtdxIxmVexIjMDDw6aefnjRpUn5+/v/7f/+vc+fOR44ccWVmqN4Tf/lJ/O1Xce9P4slkd+eCkBI5Pgd96NChjRs3Hjx4cPz48Zs2bYqKisrIyBg/fnxqaqqL80MN1vf38i/xVbtyIKXISAHuWHnbYwZAajleUAgsk3v5qt1mpfe9bbUWCqJdu0OU0pSiImDYskEoQEyAfz8/3ypljpBDjgv0kiVLZs2a9emnn9q+9W7duvWECRNcmBiq58SEfdIDmn9P/N8+hmW5IcOlFg+WBYCjRtNRY3UueDpZZCyv5ZSxuILNAOC2xeqw3fn3AoBBqWe9OG5928gJjYOdD4VQWY4LdERExLPPPmt7OmjQoN9++23x4sWuygrVZ/TWzTJNVPj1Z+HgfvVTz7Jdui8PD23vqSux8mq12lEAezyhxYQAQIrRKBBo76UTKE0zFQMwAGClpJu31/kiI8swnX28AeBicTEDTLSPNwBoAIoJLaHknKn4rtXaWK3u4u2lYxlfjlM7uopJei9K6ckiIzB/BbECAMCv+QUZJSXdvH2+79iukRqvYUU1Zf8Z2rZt2+rVq9PT00+fPi21WCwWvN4OyUj9/Czh9CkoKBBPJoPVAhwHPM8Gh6jipkgXeLT11L0THlqNK1LfuZEpUCpdByI9VjGM1PKvy1coy74T2dpus9L7Hiww9PfzdeYyElEU/+/6TeA4u42fOH8xo6REyzJhHvJcNYwaOPsCHRcXN2bMmNmzZy9ZskRqYVlWmjUJIXn4Bar6DwEAcuE8FQRGraU8D74+TPMW7s4MIWWxL9ATJ0785ptv8vPzFyxYULp9z549LswKoepYGtrS4WMAeK1JY9sNXXYv2RodtpfnjRbNyt5cNyY4cHuuPJcMIgRlC/Q///lPAJg3b547kkENi2bhmwAgnjgqfL/N3bkgpET2BXrlypUrV64su11MTIxL8kEIIXSffYHGY2eEEFII+wK9YcOGb775ZtSoUXbteASNaheu9o1QGXgOGrkdAwDk5nXxZDL30CPuTgYhBbEv0AMGDACAhx9++Kuvvjp//jzLsl26dJkyZYo7ckMNAhvWmvH2ocYiKMhzdy4IKYvjO1AmTJhw8uTJ6Ojo6OjoEydOPP300y5OCzUcTHAjpg0ueomQA47vRi0sLNy/f7/0eMaMGUOHDnVhSgghhADKO4Ju3bp1Ts79aW0LCgpatMBbvBBCyNXsj6Cl6zdMJlNERMRDDz3EcVxKSkp0dLQ7ckOozrtuthBKW+vkWQsRNTROXQednZ3tkmQQqm+2ZOeUnZgJISfZF2jb9c6pqal6vR4ALBbLihUrJk2a5OrUEKqbBErzH1zWXUKIWOqp4+0FQcAldJEjjr8knDt37uHDh7Ozs8PDw9PS0pYuXeritFADRG5eZ4sKGZ86vxbJscKiwMRjpVtW3MyqeBcWYH27yGlNQ2ozL1T3OP6SMDk5OTU1dcCAAUlJSQkJCXfv3nVxWqhBYVgWAMjFNPHAL+7OpUY6eHpqWKYaOxKAI4VFsueD6jrHR9AqlQoALBYLISQ6OvrNN990bVaoYeH6DiB3sujtLGq1ujuXGnnIxzu/by8LJSUiyRUEAPj89l2R0tnNmwJAsIrz5Liye7197caHWXgMhBxwXKCjoqJWr17drl270aNHR0RE4JeEqFYxzVtynboKtys5D1AneHKsJ7ABKmim1QBAE41aoLSrl2cFu+hwxSJUDscFet26dXq9PjAwcOPGjbm5uTt27HBxWgghhBwXaKPR+PXXX9vm4ggKCnJxWgjVD1VapQUhO44L9IQJExo3btyvXz+WZZOSkvbt27d7924XZ4YaIquVXLtKiwqBAgQFU6sFKGUC7Y8PaH6e8+1SI6jUFeyIkDLhXBxIQWjePfHALzQ/DyiwXboTrZYSoho6wm4zknLc+XapEXr2rWBHhJQJ5+JACCGFwrk4kJLc09MCDnieAjBnUkSWBQr8jQy7rWjePefbpUb11csMy4oF+aU3YJq1VI0YDUx1rlxGyAWcmotDFEWXJFNPJCcnb968uUmTJgAwc+ZMLy+v+Ph4hmE0Gs2cOXM4R1fCNnAk7azwv/0AQK0WWyOTnyfd/kzz9A73qlI7m6cHAPuA6Ze4fgMYX7/qJo5Q7XI8F4fFYvnvf/978uRJhmF69+6NE3FUiclkGjNmjG1Wk127dvXs2XPgwIE//PBDcnJy79693ZueArFt2mnmLBBTkmn6JVpYRHkzUIb6+rFqNag1bLuOoGIZjQcwlFosAECvpgOlTGRbACi3nVNRUbjfSERr46Ycx3J3bgPDMJFtgVJx/x4AAELc2HGEKub4S8LZs2eXlJQMHDgQAHbt2nXixImPPvrIpXnVZSaT6cSJE6dOnfLx8Zk2bdrNmzfHjh0LAJGRkRcuXJAKdGJiYlpaGgAYjUaz2exkZFEUKaWMTP8kJ4TI+G8jQgil1Pm+2AsIgiEjGJUGMq5AQT4AiB06M2o1JYQ8OsBuW8ZqpYRA70edaZcarb36qVQqJunw/Q1Ekdm/BwAsFgtUPWdCiNTfqu7okPRTEEWx+qP3d4IgyBUKAAghPM/L1VnygCzRpN8IGTtbG0PnZGcdDovjAp2RkXHo0CHp8ezZswcPHlyTLBua3r179+rVKzg4eMeOHQkJCbZ2Qoittl68eFF6qUOHDjzPOxlZ+vk5v33FKKVy/Z5I0aDGuXGiyFLKAKUAlFJRFIFSsUxMrirtUqMU7a+AoqgBuJ9w1XOmD1R1R4ekn4L0yyxXQLlCwYM/5HJ1Vvaho5TK21kZo1FKBUFwvkCXHRbHBZrneUEQpBk57n+skdOysrKkW3t0Oh0hJDw8PD09vWXLlunp6eHh4dI206dPnz59OgDEx8f7+Pg4GbmoqEilUul0OlnyFEWR53kPD3nmkhdF0WAwON8Xxx77h+2hyWTy9PKqdLPK2x/7BwAIhYUajcbDtoEoSme7vb29marnLIqiIAharbaqOzqkVusBQK1W13T0HjCZTF7lDV3VFRYWarVauTorFSyNRiNLNJ7njUajXOMGcg+dwWDQ6XROdpYQUvYLKscFeujQoX379pUOnH/77bfHHnushok2KL6+vmvXrvXx8WEYRvrS9aOPPjp27JiXl9eTTz7p7uwQQnWG4wK9bNmygQMHHj16lGGY//znP/3793dxWnVaWFjY8uXLS7csXLjQXckghOouBwWaELJhw4aZM2cOGjTI9Qkh1MDJu4xh2WgVx6/Su9cwVaUt2HjdbDGaLZ3kO8VRcw7uJGRZdufOnXl5ea7PBiG0JTvnv9k5tRet4vhVevcapipvT2tuS3bO13n57s7ibxyf4mAYJjIyskOHDrZvkEpfjYAQkl2WxZqQXwAAV0vMIlDpcfWYzWYP6/1LEcpGqzh+2VeLi4vVZotarXZm40pJF4RIFyDUvKfSVXHeIM+Fp7k87y1LIPk4LtCvvfaaw58HQqiW7MvL31fq8G3z3VwZg5eNVnH8Kr17DVOVt6c18YiPd4yvskq04wL90EMPffnll7b5oKdMmeLitBBqOB4LDNieoy+kkMNbAYBSAGAYhgJAI7Xaq+pzA1BKSwi1i0YpAANMOfGLRVLeu0vX75e+PaqCjZ3JrVgkuQIvV08ppWzNlqSxdSfFaDptKv44OxcA1reNnNA4uCZhZeG4QD/11FNhYWE9e/YURRHng0aoVvXy8T4aFan19cvheQBYm3VHoHRei2YA0EitqkbZMplMjIfOLlqJSChQT45zGL9YJOW9e9nroCvYuNLcBEEw8kIBw8jSU+k66ICAgKruWFrp7pRYra+Fh0r51CSmXBwnodFoPv30U+nxrFmzpHu+EUK1x5NjwzgtAPirOIHSMI8a3RhSQTSH8av07jVM1ZNj/TWaau8uu9Ld0RLO7fmU5vifBt7e3iaTSXqcn5/ftGlTF6aEEEIIoLwj6Nu3b7dq1ap79+6EkJSUlK5du0rzRO/Zs8e16SHU4Mi7jGHZaBXHr9K71zBVpS3YuDS0pe3AVCEcF+hFixbhtMUIIeRejgv08OHDXZwHQgghOzW6PAUhhFDtwQKNEEIK5XjR2LLw60FUL4lHExmPv8+vrdZwDz0MHvJMuo1QTTi1aKxe73h1ToTqKoYBhgFKxd8OlH2RGgtVwx0fqSDkSo4XjQWA1NRUqS5bLJb33ntv4sSJrk4NodrDsqrHxpDMG3bN9PYtqs8Fq9UtSSFkx/FVHHPnzj18+HB2dnZ4eHhaWtrSpUtdnBZCtY3rP7jslaTC7h1i4u+uTwYhhxx/SZicnJyamjpgwICkpKSEhIS7d++6OC2EEEKOC7Q0W6vFYiGEREdHp6WluTYrhBBC5ZziiIqKWr16dbt27UaPHh0REZGdne3itBBCCDku0OvWrdPr9YGBgRs3bszNzd2xY4eL00IIIeS4QP/666/Dhg3jOG7GjBkuTgghhJDE8Tno1atXt2rV6pVXXjl58qSLE0IIISRxXKATEhJOnTrVpk2b+fPnd+7cefny5S5OCyGEULlzcTRu3PiFF15YuXJl79693333XVfmhBBCCMor0ElJSfPnz2/duvWSJUv69Olz+/ZtF6eFEELI8ZeEL7/88qRJk5KSkpo0aeLihBBCCEnsC/TEiRO/+eabkJCQhISEhIQEWzvOZocQQi5mX6D/+c9/Qjlz2iGEEHIl+wI9YMAAAOB5XroO2h0pIYQQAsDroBFCSLHwOmiEEFIox1dxwIProKOjozds2PDuu+++8cYbrkyr4RAEwWQyOb8xIYQQIstbS6FEUZQlGqWUUup8XyrF87yM0QRBoJRW2lmW51kAnuctDt+6IB8ohYBASikhRBCE0o3Vzo0QouShE0XRmpMtMAwEBMrSWSlDWXKTPsNKHjqz2exkZx3+Mjou0ElJST/88MOPP/4YERExadKkVatW1TRTVA6VSuXl5eXkxoQQlUql08mzXJ4oijzPe3h4yBXNarU635dKmUwmGaOJoqjRaCrtrKBWiwBqtVrl6K3Fo4mUENXQEaIoCoKg1WpLN1Y7N0EQBEFQ9NClneVYVjV0hCydJYRoNBpZcuN5XslDJwiCh4eHk50lhJT92g+vg0bIHi0qpFmZDtoLDUApzcoEQkAQqEbzt8ZqE0XWZKLFxupH+DumpITK9FccAFiTCYoKgWVpVmZVO8s0agwarVyZNEAOCjQhZMaMGTNnznR9NggpATmdYj2dUt6r4rE/pQdWR43V4/H3aDXEyRpNDQAAYqk+Ot9ZpnGI5tUl8uXS4Dgo0CzL7ty584knnggMrP6ZJoTqIjaiDTl1gtqd5RfFUsvIUgAGAEClAukcdOlGjQaqdXEqpZRhmOrlXIvR/uo4/fsLTnSWEjCbqT5XhjQaMMenOBiGiYyM7NChg+2cXem7ChGqr9gOnTVvdrZv5a3UaAQAMekwEML1GyiKoshx0plFWyMAMF7eUPWzq4IgFBUVBQQE1DT7B2Q7kcpbqdFoNBo1J49ygsBG9yInk4HjuAFDoLLOUn2OdeU7MuTQsDku0K+99pparXZxKggplFrDBAQCAOOho4QwAYGMKIIgMFpt6UZ3Z1kL1BomIJByKtB5AcuyYa3p1fR621lFclygs7KySj8lhPTr188l+SCEELrPcYH+/fffpQeCIBw/frxjx46TJ092XVIIKRI3ZLiTjfXNwBiVVgsNpLNK4rhAb9iwwfZYFMW5c+e6Kh+EEEL3lbuiig3Hcbdu3XJBKgghhEpzfAQ9duxY2+OcnJygoCBX5YMQQug+xwV69uzZ0gOGYQICAqKjo12YEkIIIYCypziKioo+++yz2NjY2NhYvV6/bNmyjz/+uKCgwC3JIYRQQ2ZfoGfPnn327FkAyMjIeOWVVxYvXhwREfHKK6+4IzeEEGrQ7E9xnD179syZMwCwY8eOSZMmjRo1auTIkd26dXNHbggh1KDZH0F7e3tLDw4dOjR48GAAYBhGxlkCEEIIOcm+QIuieOfOnczMzMTExIEDBwKAXq+nlDrYFSGEUG2yP8WxaNGiTp06WSyWd955x9fXNz8/f9CgQbNmzXJLcggh1JDZF+hx48bFxsZaLBZ/f38A8PX1ff/990eOHOmO3BBCqEFzcB20TqezLarEcRxWZ4QQcovKb/VGCCHkFligEUJIobBAI4SQQmGBRgghhcICjRBCCoUFGiGEFAoLNEIIKRQWaIQQUigs0AghpFBYoBFCSKGwQCOEkEJhgUYIIYXCAo0QQgrleFVvVBMGg2HNmjUajcZqtS5YsODcuXObN29u0qQJAMycObNRo0buThAhVDdggZbf9evXR40a1aNHj61bt547d85kMo0ZMyYmJsbdeSGE6hgs0PLr2rUrAAiCcPXq1ZEjRyYmJp44ceLUqVM+Pj7Tpk1Tq9UAkJKScvXqVQAoLi62WCxORhZFEQCc375ihBBRFGWMRimVKxoAyJgbABBCBEGQd+hkCQUAoigqfOh4nq/yblar9F+7TKTOyrWKXp0YOic7K/0G2TViga4VBoNh3bp1cXFx/v7+vXv37tWrV3Bw8I4dOxISEkaMGAEAiYmJu3fvBoDevXubzWYnw5IHZElS+j2Rq9BIny3n+1Ip6XdP3mgydlb2aIodOqnKVLWzjMWiBYAyHwmps4IgyJIbpZQQovChc7KzDv/qY4GWX1FRUXx8/OzZs4ODgwEgKysrKCgIAHQ6na22zp07d+7cuQAQHx/v5+fnfGSVSmVb76aGRFHked7Dw0OuaAaDwfm+VMpkMnl5eckVrbCwUKPRyNhZQRC0Wq0s0QRBKCoqUvLQabXaqnaW8hbpENquX4IgEEI0Go0sufE8bzQaFTt0BoNBp9M52VlCiEplX5CxQMtv9+7dubm569atA4CYmJiQkJC1a9f6+PgwDDNv3jx3Z4cQqjOwQMsvLi4uLi6udMvy5cvdlQxCqO7C66ARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaIYQUCgs0QggpFBZohBBSKCzQCCGkUFigEUJIobBAI4SQQmGBRgghhcICjRBCCoUFGiGEFAoLNEIIKRQWaIQQUigs0AghpFBYoBFCSKGwQCOEkEJhgUYIIYXCAo0QQgqFBRohhBQKCzRCCCkUFmiEEFIoLNAIIaRQWKARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaIYQUCgs0QggpFBZohBBSKCzQCCGkUCp3J9DQCYJgMpmc35gQQgiR5a2lUKIoyhKNUkopdb4vleJ5XsZogiBQSmXsLCFEEARZohFClDx0oihaLJaqdpYpKeEAAMAuE+nTy/O8LLlJn2ElD53ZbHaysw5/GbFAu5lKpfLy8nJyY0KISqXS6XSyvLUoijzPe3h4yBXNarU635dKmUwmGaOJoqjRaGTsrCAIWq1WlmiCIAiCoOSh02q1Ve0sLTFZAQDALhPpIEOj0ciSG8/zSh46QRA8PDyc7CwhhOM4u0Y8xYEQQgqFBRohhBQKCzRCCCkUFmiEEFIoLNAIIaRQWKARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaIYQUCgs0QggpFBZohBBSKCzQCCGkUFigEUJIobBAI4SQQmGBRgghhcICjRBCCoUFGiGEFAoLNEIIKRQWaIQQUigs0AghpFBYoBFCSKGwQCOEkEJhgUYIIYXCAo0QQgqFBRohhBQKCzRCCCkUFmiEEFIoLNAIIaRQWKARQkihsEAjhJBCYYFGCCGFUrk7gfqvuLg4Pj6eYRiNRjNnzhyO49ydEUKobsAj6Fp34MCBnj17vv766y1atEhOTnZ3OgihOgOPoGvdzZs3x44dCwCRkZEXLlzo3bs3AFy8ePHWrVsAYDabeZ53MhQhRBRF57d3cTQAkCuaFFDeaIodOlEUKaWKHTpKaXU6KwjSf+12lLezgiDUm6EjhFBK7RqxQLsOIYRhGOnxjz/++MMPPwDA8OHDTSaTkxFEURRFUXjw0a8hSqm8H25KqfN9qZRUBOWKJoqijL970i+S1WqVK5ryh66qnWVUambiFGAY8vd+SUNn+0WoIUopIUTJQ2c2m53srMO3xgJd68LDw9PT01u2bJmenh4eHi41Llq0aNGiRQAQHx/v7+/vZKiioiKVSqXT6WRJTPrb7uHhIVc0g8HgfF8qZTKZvLy85IpWWFio0Whk7KwgCFqtVpZogiAUFRUpeei0Wm11OhvcqGybIAiEEI1GI0NmADzPG41GxQ6dwWDQ6XROdpYQolLZF2Q8B13rhgwZcvz48eXLl2dnZz/88MPuTgchVGfgEXSt0+l0CxcudHcWCKG6B4+gEUJIobBAI4SQQmGBRgghhcICjRBCCoUFGiGEFAoLNEIIKRQWaIQQUigs0AghpFBYoBFCSKGwQCOEkELhrd5uptfr4+Pjndw4Pz9fo9HINZmLNBOYXAsICIKg1+ubNGkiSzQpYNm5Y6rt3r17Op3O09NTlmjyDh3P83l5eSEhIbJEA2UPnTSpplxDZ7FYDAZD48aNZYkGcg+dXq/39vZ2coouSmlAQICDVlRXPPPMM2vWrHF3Fo6lpaWFhIS4O4tyjR8//vPPP3d3Fo6dOHEiPDzc3VmUKzY2duvWre7OwrHDhw937NjR3VmUq3///jt37qxJBDzFgRBCCsW99dZb7s4BOaugoKBt27a2SaUVRRCE4uLiYcOGuTsRx/Lz8zt06NCqVSt3J+KAIAgWi2XIkCHuTsSxvLy8Ll26NG/e3N2JOMDzvCiKAwcOdHcijt27d69Hjx41Oe/H0DKLrCCEEFICPMWBEEIKhVdx1A3FxcXx8fEMw2g0mjlz5sj1JXi1JScnb968Wfq328yZM728vEqnZ7FY3JItIeTHH3/cuXPn5s2bocyg2WXl4iTtclPaABoMhjVr1mg0GqvVumDBAgBQztDZ5Xbu3DlFDd2tW7fWrVvn6elZUlLy6quvqtVqGYcOz0HXDb/88ktERMSkSZNu3bpVVFTUsmVL9+Zz+fLltm3bTp06tX///l5eXnbpnT171i3ZFhQUaLXajIyMoUOHQplBs8vKxUna5aa0Abx06VKbNm2efPLJW7duiaJ4+vRp5QydXW4mk0lRQ5eXl9e/f//hw4dfvXrVy8srJSVFxqHDUxx1w82bNyMiIgAgMjLy+vXr7k4HTCZTYmLiypUrP/vsM57n7dJzV7YBAQEdO3a0Pa04KxcnaZeb0gawa9euPXr0EATh6tWrbdu2zRK6wQAAIABJREFUVdTQ2eWmtKELCwuzWCyLFi3Kzs7u0qWLvEOHpzjqGEKIXEvW10Tv3r179eoVHBy8Y8eOhIQEW7tdegrJtuKs3JKkAgfQYDCsW7cuLi6u9CLZChm60rkpcOiaN2/+/vvvb9u27dChQ07m42R6eARdN4SHh6enpwNAenq6Ei6zy8rKslgsAKDT6QghdukpJNuKs3JvkkobwKKiovj4+Oeffz4yMhIUNnR2uSlt6H744YezZ88CQEBAgMlkknfo8DK7uqGkpOSjjz4ihHh5eb300kss6+a/rNevX1+3bp2Pjw/DMPPmzQOA0ulZLBa3ZHvx4sXt27enpaV16NBh8ODB3bp1qyArFydpl1vTpk0VNYDbtm07duyYdM90TExMly5dlDN0drmFhIQoauhycnI++eQTDw8Pq9X66quvsiwr49BhgUYIIYXCUxwIIaRQWKARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaoSqYNGlSZGRkaGioWq2OjIyMjIxct27doUOHJk6cWO2YhJDRo0cnJSWVbpw9e/bGjRurF3D+/Pnr16+vdj5IOfBOQoSqYOvWrQBw8eLF2NjYK1eu2NoHDBhQ7ZiffPJJ+/bt+/TpI0N+AACwYsWKbt26xcbGKnP+a+Q8PIJGqKYSEhJGjRp16NCh4cOHz507d/DgwdOnT1+5cuWUKVN69uyp1+sBYPny5YMGDRo8ePC0adN4ni+9+6pVq+bPnw8AoihOnjy5b9++Tz/99O3bt6WWGTNm9O3bt3fv3nPnzgWAvn372u4njo2N3bt3b05OzrBhwwYOHNi1a9fPPvsMANRq9ezZs9esWePicUCywwKNkDw4jjtz5szKlSsPHDjw3Xffde7cedOmTd26dfv5558TExMTExMPHjx48OBBjUazadMm214XL1709/eXJs/cu3dvVlbWn3/++dVXX50/fx4A9Hp9VFTUn3/+eeTIkcOHD6emptpOfeTl5V2+fHnYsGHbt2/v2LHj77//npycLN0DDQDDhw//9ddf3TAKSFZ4igMh2XTq1Emr1QJAYGBg9+7dAaBRo0YFBQWJiYnnz58fNGgQABiNxkaNGtl2yczMtJ2IOH36dK9evQDAw8Ojb9++ABAUFHTv3r1x48Zptdrs7Gy9Xv/kk08uXbrUZDLt2LFj0qRJHMfFxMSsWrXKYDDExsbOnj1bChUaGpqZmena3iP5YYFGSDal51+3PaaUMgzzxBNPrFq1quwu0qu2x7Z2URQBYOPGjampqbt37+Y4TjpJ7eHhMX78+B07dnz//ffSCY02bdqkpaUdOnRox44dK1asOHHihBJmEESywFMcCNW6fv367d2712w2A8CqVatSUlJsL7Vq1ermzZvS444dOx49ehQATCbTH3/8AQA5OTkREREcx506derChQtWqxUAZs2a9fHHHzMMI02Htn379tTU1GHDhn388cc3b940Go0AkJmZ2aJFC5d3FMkMj6ARqnV9+/Z99tlnBwwYoNVqW7Vq9eKLL9peioqKKigoyM3NbdSo0ejRo7/99tuePXu2aNGiX79+hJC4uLjx48fHxMR06tRp8eLFixYteuSRR6KiogghU6dOlSK0b99+5syZKpXKarW+9tprPj4+AHDgwAHFrrCOnIez2SHkZmvXrr1169b777/v5PYZGRmPPfbY2bNn1Wq1ww0EQejevfuePXtCQ0PlSxO5AZ7iQMjNXnjhhfPnz9vdqFKed999d/z48V9++WV51RkAFi1aNGfOHKzO9QAeQSOEkELhETRCCCkUFmiEEFIoLNAIIaRQWKARQkihsEAjhJBCYYFGCCGFwgKNEEIKhQUaIYQUCgs0QggpFBZohBBSKCzQ7hccHNytWzdBEGwtEydOdDh3sJNu3brFMIw07aS8vv/++5ycnGq/de0lhlC9hAVaEQoLC1evXu3GBCilhJBKN3vzzTcrKNCU0iZNmty5c8fLy0vW7BBqoLBAK8LKlSvfe++9q1ev2rVfuXKFYRhponcAiImJWbt27Z07d1iW3bhxY0xMTIsWLd5+++3169fHxMSEh4evXLnStu/XX3/drl27Ro0aTZ48uaSkBACSk5P79esXFRXVqVOnTz/9FAAyMzM5jvv444+Dg4MvXrxY+q0TEhK6d+8eFhYWHh7+4YcfAsATTzxx8eLFf/zjH1999VXpLUsHOXjwYNOmTU0mEyHklVdeiYiIaNOmTbdu3fbv32/Xta+//josLAyXZUKoIhS5W1BQ0Llz55YsWTJkyBCp5amnnlq5ciWlND09HQBKSkqk9iFDhsTHx+fm5gLABx98QCk9cuQIx3ErVqyglJ46dUqj0ZjNZqnqzZo1SxCEvLy8yMjINWvWGAyG4ODgrVu3UkqzsrKaNGny22+/SYfDb7zxhiiKpVOyWCz+/v4//PADpfTMmTMcx6WlpVFKAeDs2bN2+ZcOIr11UVHR/v37w8LCiouLKaWJiYmTJ0+mlNpePXDgQIsWLS5dulSLw4pQ3YdH0IpAKV26dOnNmzdLr/dcsccffxwA2rdvL4ri+PHjpcdWq9V2CmLu3LkcxwUEBDz++OO///77b7/9ptPp4uLiAKBZs2ZxcXHffPONtHjdpEmTWPZvnwSNRnP9+vVx48YBQOfOnZs1ayb9qXDIYZBmzZrp9fpNmzbdvn27b9++pft18uTJ559/fs+ePW3btnWyswg1TFiglcLDw2P9+vULFiyQDpArJa1sJK1M6u3tbXssLTYKAE2bNpUeBAUF5eXl3b17Nzs7O+yBbdu25eXlSRs0bty4bPxNmzb169evV69evXr1ysnJqfQMtV2QTp067d69++DBg507d+7evXvpUxyTJk3ieT4kJMSZbiLUkOGahAoycODA0aNHz5s3z7ZYhlRz6YNFFQoLC52PptfrAwICACAvLy84OLhZs2atWrX6/+3deVgT9/Yw8JOEBAKGHQRFFkEMq7JUrKgsbmjlutxab0Gr1raXW5cfrdentre37e/26as+xW7YFqlVLFfrUsG2VnxFqXhRkbqAbBFQLKAVhJgQAmSZ+b5/jKa5YSkmo0zfns/DH5khc3JyvnAYJjPfMdkRbm9vhwe7wMYKCgrefvvtn376KSAgAABGjx79my/XN0h8fHx8fDxN0/v27Vu8eDHzWgBw5MiRnJyc1NTUwsJCkz13hJAx/PXgloyMjJMnT5aVlTGLHh4eQqGwqqoKAMrLyysrK4ceKjs7GwBUKlV+fv7MmTPj4+OVSmV+fj4AaDSa9evX//jjjwNte/fuXXd3dz8/P0LIhx9+2NXVxZwbJxQKDfvdg9u9e/cLL7yg0+n4fH5YWJjxSYRSqXTbtm1yufxf//rX0N8OQn9A2KC5xdnZ+aOPPqqrq2MWxWJxRkbGsmXLZsyYsXPnzrlz5xp3uoHodDoejxcSEjJhwoSgoKCpU6euWrVKIpEcPXp027ZtAQEBUqmUoqgnn3xyoAh//vOf/f39AwICJkyY4OzsvHbt2ldeeeXMmTMpKSnz588fSmN9+umnNRqNv7//2LFjn3vuuf3794vFYsN3ra2tDxw48OGHHxYVFQ2hKgj9QeE9CRFCiKNwDxohhDgKGzRCCHEUNmiEEOIobNAIIcRR2KARQoijsEEjhBBHYYNGCCGOwku9h5larTbMnvGbmJPW+15UbTZCyB8nGrBaOnZxvHScrRtw+80+bDRra2tra2vjNXihyjDbuHFjaGjoEJ+s0Wj4fL5hpg4L0TRN07SVFTt/pGma7u3ttbW1ZSUaAGi1WpFIxFY0jUYjEAhYfLPslk6j0RhfaWkhjpeOEMJMMmM5iqK0Wi1nS9fb2ysUCof4Zmma7unpefnll41X4h70MPP29l6xYsUQn6xSqaysrNj6caQoSqfT2djYsBVNqVQ6OzuzEg0A1Go1i3dm6ezsFIlELL5ZvV5vsrNjNr1er1KpmJmtWMF66fru2ZlNr9fTNM1WE9TpdF1dXZwtnVKpFIvFQ3yzNE1nZWWZrMRj0AghxFHYoBFCiKOwQSOEEEdhg0YIIY7CBo0QQhyFDRohhDgKGzRCCHEUNmiEEOIobNAIIcRR2KARQoijsEEjhBBHYYNmH33zhv7/HjVZqT9yiPxy+2FDnThxIjIyUiqVhoaG7tu3j1l56tSpyMhIJyensLCws2fPspCxWY4dOxYeHu7o6Dh+/Pi9e/dyKjdjTA2jo6OjoqIMNeSOY8eORURE+Pr6GpeRI5ghHjNmTFhYGNdyA4DS0tInnnjC29t7/PjxBQUFw53OfyktLY2IiPDy8goNDbUoN2IuSlaj2fKOdlcW80X93EjJajT/eoNWdxFCaMU97dd7iE6r2/+V9t+7tF/t1O3/iuj1hBD9xQva7Eztzs+0n31I1VQSQkhPtzb3S+2/d+kO5BKKsmiRovSnCzX/+/r9LGlad+SQNvdL7a7PSU8Ps05/uUy7fVvf19X/VKrNydbmZGs+2qI/d4Zub9Pu+lx3cK82J5vuUj1420bxTV7rAe0X24lWc//p5Zfur1Wrtbt39C3jJ598MlCF7927Z2dnV1paSgiRyWRisfjKlSuNjY0SieTIkSOEkL1797q6uqpUqoEiDE6v1/c8qMnDam5ulkgkRUVFhJBTp06JRKKGhoba2lq2ciOEdHV1mb2tgaGGSqWyoqJCLBbX1dVZHlav1/f29loehyljYWGhXC5nytjY2Gh5WFZKZxhipVJ5/PhxtnLT6XQajcbyON3d3SNHjvziiy/a29uPHz9ub2/f2tpqeVhWSsfktmvXLrlcfvTo0SHmRlHUp59+arLSoj1ofugE4aq/Ml98b18A4PmPo3741vAE+udGsLUTpq4SLl/NDwol3Wq66We6/KLw+b8JV/9NuPpvdHMTAFA/lQpCwoWpq3juI+maKksWibqL5+3Hc3O/n0DDNRCLhcuet5r7J0JTAEC6VKT+GohEfV9XEB0jXPGicPlqnoubIDqGbqjnh06wWpLC8xxNmpvu/z0zim/yWvefcOc2z8kFhPfnr6Kv1dz/hq0tCIVEcW/o5W1sbBQIBDExMQAwfvx4T0/Purq6c+fO+fn5LViwAABSUlLc3d0LCwsfduAsx+Pxvvrqq4SEBABITEx0c3O7fv16aWkpF3IzZlzDwMBAT09PmUw2vCkZY8oYHx8PD8rY0NAw3EndZzzE8fHxnMoNAIqKihwdHVesWMHn8+fMmTNp0qTvvvtuuJO6j8lt1apVfD5/1qxZluQ2pOlG29vbXV1d+66nqyt17XeZx8JnVwAA38ub3JPT9dd47iMBgO/lQ50p0n93mD82gB8oBRsxXXZOEBUDzASpImur2fMAgLTe4Y8PBgDe6DHk50aiUJi9yA8N50nsDRPgk+Ym0Gr1efuBL7CavxAAqILvBXOS9Qdz+74uhIYDAHXhrGBiNAhFgogobfZ2+lotaDT8mUlMQJ7E3hDf+PGvNbnVwhvtBQDUxVK6Tkaabur25fBHewniZvJHe5HbLTxHJwBQKBRqtRoAKIoaaML+wMBAiURy9OjRuXPnXrp0SaFQREVFnT9/3ngTW1tbmUw29Cn/jVEURdO0edt6eHgkJyczs/EeOHAAACIiIm7dusVWbgBgdm7GDDWcNm3axYsXFQpFdHS05WEtKZ0x4zLm5OQAQGRkpOVh2c1No9EcPHiQxdxYSa+2tlYqlTKzS1MUFRgYWFNTw5HSMblRFEUIoWl6iLnRNN135YANury8vL29HQA0Gs3rr79+9erVvs/hh4RZPbXQNGJSsm53ltXTzwIAWFsLV6URxT26sUG3e4fVgiWEpnmEAABdL6PO/EjUXcK/PPfrxoQA8NhZZNZRep6TsyA2jvrPj9Tln3g2Yp63L8/BwfSdGLYlhC6/JEz7HwCgSoqtpsbzJ0ZRZ4qoKxcFUZMGqtV/6enmSewBQBA9WRA9WX/w31bPLLv/LRsx6elhHubk5DB/VKdNm6ZQKAYK9sEHHyxdutTW1razszMjI8PBwSEqKurmzZu7d+9OTk4+duxYTU1NYmLiIBEGRwjpeZCSGXJzczds2ODm5paVlcXn8ydPnsxublqt1uzcDExqKBKJzE7JRHd3NytxjMtICLE8PbZK9yhyY0tHRwefz2eOSCgUCh6PJ5fLOVI6JjeFQkEIUavVQ8yNpmmdTmeysv8GvX79+jNnzrS2tvr5+dXU1Lz55psPkZ21tWB6IlV0AgDoa7UgFPLHBgginoCuLrqpUSAN0f9whB8azh8n5Y+T6vfngl7P8xxNN/8scB9Jmn/mjRoNYlvzF/8bf6QnUSoAAGxsQK+nZdUAoP/5Bmm9Q104Z/q6AHRzE2+0F/B4AEC61TyvMQAAYlvoVg/17duISe8ALa+3h2d//29Denp6eno6AGRmZrq4uPT79Pr6+rVr1xYXF0dHR9+4cSMhIcHPz2/GjBl5eXmvvfbaW2+9tXDhwri4uFGjRg0UYXCWT9ifnp6+Zs2a8+fPL1269NChQ8HBwWzlBixNnW6oYWBg4K1bt5KSkkJCQhITEy0My+6E/WvXrk1JSZHJZEuXLj18+PCUKVMsDMjirPPp6enLly8vLy9ftmwZK7mxNWG/u7t7RUWFvb09M2G/RqNxc3Mz+4fNgJXSMbm5uLgwE/YPMTeapvveLKn/Y9BlZWXl5eVxcXHnzp07efLknTt3+o9YVaHbvYP5oisu/xo0OIxotQDAG+1FnS3W5e7U5X5JtzQJJkbxxvjwY6bovvxct3uHbs8XPK8xPM9RgugYWlaty/2SyOV8aYgli/TPjbo9X9Btd3R7vqCrKvgh4eSXW7p9OXR1pSAi2uqZZcwXb6SHIGaKybYAQFpv85zvH8wRTJlOnS/Rf/M1fa1GEBVDWpp0+78yjq//4Yjxa91/76O8yO1bhlL8uvvMHP0Y5TX4IBkrLi6WSqXR0dEAMHbs2NjY2OLiYgCYPXv2lStXbt++vX379qqqqsjIyKHHZEt1dTXzH4BQKJw+fXpcXBxzuJkLuRkzrqGfn19sbOypU6eGNyVjfct44sSJ4U7qPuPcpk6dyqncACA4OLimpsawWFlZGRISMoz5GGMzt34/T4yNjSWELFy4kDmMMnfuXDM/zvxD0mZnEp3WdO3Dn8VRWlrq4ODAnHXQ1tbm7e29Z8+eO3fu+Pr6XrlyRa/Xb968OTQ0lBkjM1hyFkdJSYlEIrlw4QIhRCaTeXh45OfnNzY2spUbYenzdEMNlUplc3Ozt7f3wYMHLQ/L1lkcTBnPnTsnl8uZMn777beWh2WldIYhViqVV69eZSs3ts7i6O3t9fLy2r59e1tb2zfffOPk5CSXyy0Py0rpmNw+//zz9vb2/fv3DzG3fs/i6L9Br169OiMj47XXXps3b966desiIyMtT/qPg7rRoDv+vclKXf5B+nZL3ycP0qAJIdnZ2UFBQf7+/uPGjXv33Xc7Ozu7u7tzcnJ8fHwcHBzi4+MbGhrMztOSBk0I2blzZ2BgoIODg6+v79atW/V6fUdHB1u5EZZ+VciDGvr5+QUEBLz77rusxGSrQRNCdu7cOW7cOHt7e6aMrMRkq3TMENvb2/v4+LCVG1sNmhBy6dKliIgIsVgcHBx8+vRpVmKyVbpLly5FRkaKxeKgoKAh5tZvg+7/rt4URbW3tzs7O+fk5Ny9ezc1NdXHx8fMXXQ0qMzMzHXr1g3xyXjTWLPhTWPNhjeNNZsZN40d0l29n3/++T179gDAiy++CAAJCQk//vijxdkihBB6CKYNet++fdu2bauvr6+ouP+pl0aj4fPxinCEEHrcTBt0SkrKggUL0tLS/vGPfzBr+Hy+ra3tY08MIYT+6PrZNbazs8vNze3t7W1paWlpaamvr583b97jzwwhhP7gHsGFKgghhNhg0YUqCCGEHp3+G7SVlRUAaDQamqajo6ONr4pBCCH0ePR/iEMqlW7btm38+PHJycn+/v6tra2POS2EEEL9N+gdO3YYX6iSl5f3mNNCCCHUT4OmafrLL7986aWX4MGFKgghhB6/fo5B8/n8/Px8uVz++LNBCCFk0P8hDh6PFxAQEBwcbJi74OTJk48xK4QQQgM06I0bN/adOhohhNDj1H+DzsnJYSZLYuBkSQgh9PjhZEkIIcRRQ5osydvb+7EnhhBCf3Smu8YqlSo3Nzc3N1cqlV68ePH5559/77332Lp7MUIIoaEzbdBpaWmVlZUAcOPGjVdeeeWNN97w9/d/5ZVXhiM3hBD6QzM9xFFZWXn16lUAyMvLS01NnT9//lNPPTVx4sThyA0hhP7QTPegR4wYwTwoLi5OTEwEAB6Px+PxHndev38lys5/NP5sWFxbf+Nql/phg5w4cSIyMlIqlYaGhu7bt49ZeerUqcjISCcnp7CwsLNnz7KW8UM6duxYeHi4o6Pj+PHj9+7dy6y0+W9KpXK40gOA7du3i8Xi/fv3G9ZwKj3GiRMnnnjiiUmTJhkPMUcwQzxmzJiwsDDDEHPKZ5995unpaTzEHFFaWhoREeHl5RUaGlpQUGB+IJObyE6aNOn27dtNTU2Ojo5KpZIQcvfu3fDw8IHuRFvQIR9bevGpq9XM13llZ0GH3P3shXatjhDS3KtJrbnWQ1HLa+qWVMkWV9Uur6nT0jQhJOeX1hnlVXMqqmIvX/2+vYMQotDpn66qfaZatqK2TkfTlizqaXrrzy1uJReYJGlC1tZd/3NV7byr1UqdnlmZe6dt8qWKvq+7+5fWBZU1CyprJv505dOW2/XdPfOuVq+qrVtQWdOm1TLbGsc3eS2DmeVVaj1lWOzQ6uZfrelbwEHu6n3v3j07O7vS0lJCiEwmE4vFV65caWxslEgkR44cIYTs3bvX1dVVpVINFGFwltzVu7m5WSKRFBUVEUJOnTolEokaGhpu3boFAMyPjeUsvL/y6tWrV65cGR4e/vXXXxNClEolc3EsK+mxdVdvZojPnj0rl8uZIa6rq7M8LCu3pjYMsVKpPH78uEgkamxstDwsi3f1Xr169XPPPRcSEsIMMStYKV13d/fIkSN37doll8uPHj1qb2/f2tr6m1v1e1dv0z3oTZs2hYaGBgUFvf322/b29vfu3UtISBh8Ro5Fri5Hw4KZr8n2EgBIcHTYeL3R8IRzSpWz0OpgyPjDIdL5Lk7tOt2FTtW+trvHwoKPh4ccDw8u6+wCgF13Whe6uhwIHi+1FX/fIbdksU2ne9JBEmh7/+7XJ+8pHK2svgmRbhnrqyMEAFq1upP3FHYCQd/XXenhfiQ06HCINEBss9JjZNE9xWJXl13SceF2dj91djEBjeObvBajUt3tZ2NtK/i1vM5CKzGf36zRDP1vZ2Njo0AgiImJAYDx48d7enrW1dWdO3fOz89vwYIFAJCSkuLu7l5YWDj0mGzh8XhfffVVQkICACQmJrq5uV2/fl2pVPL5fIlE8vjz6SstLW337t3Gd0DnVHoMZognTZoED4ZYJpMNd1L3GQ9xfHy8m5tbQ0PDcCf1X9LS0nbu3MnWndpZVFRU5OjouGrVKj6fP2vWrEmTJn333XfmhTI9Br1o0aKkpCSNRuPo6AgA9vb2W7ZseeqppwYJ8W17R0NPD/N4b1AgAERLRvzcqzl5TyG1tQWAJ+xHbGu59T8NN+IcHOY4OzlYCXb+0rrSw13E5wHACIHgX37eAFCt7k4aMxoAoiQjzilVzRqN2YuLXF08jW51XtbZpaaov9Zdt+Lxtvn7AsDrN26+5+ezQlbf93UXuboAwI7bd1JGutkK+Ckj3WZUVBXI73VR1Fu+Y5iAniKRIb7xY4PLqq5IyQiTlZESu8sq9RhrawBQKBRqtRoAKIqiKKrfwgYGBkokkqNHj86dO/fSpUsKhSIqKur8+fPGm9ja2spksoEiDI6iKJqmzdvWw8MjOTmZoiitVnvgwAEAiIiIaGhoEIlETz/9dGlpqYuLy+uvv/7MM8+YEZxhdm6MiIgIiqIIIUwcQohcLmcrPUtKZ4wZ4mPHjk2bNq2srEyhUERHR1selpXcDEOs0WgOHjwIAJGRkazkxkp6ABAREaHX64Gl98tgJVRtba1UKjX8+AUGBtbU1PxmWJqm+67s50pCsVhs2O8QCASDd2cAWODqkuHva7Ly/4z1mV9Z80VgAABIBIIfwoKbejVnlJ1PVdZsHzdWTwhNAAAK7ym2Nd9q1+lzgwJ/TZSA8TFvSxYZWkL72tis9/L8sOV2butdB4Fgsr1ktLVpVzVsSxP4uq39TEQYAHzc8ku616hn3d22Nd/a23p3hYf74NVgKPT6kSIRAHzUcvu4/N5Sd7dVHu4OVlYKvZ55Qk5ODvNHddq0aQqFYqA4H3zwwdKlS21tbTs7OzMyMhwcHKKiom7evLl79+7k5ORjx47V1NQkJiYOEmFwhJCeB39czZCbm7thwwY3N7esrCw+ny8SiRYvXpySkpKVlVVSUpKamurq6hoZGWl2blqt1uzcGBRFdXd3KxQK5h9GFtMDAFZOP/3ggw+effZZwxCLRCKzR9OAldIxjIeYEGJ5bo8CM8SshGKldB0dHXw+n/mpU6vVPB5PLpf/ZoY0Tet0OpOV/V/qbTmJQPCq1+j3mloA4Lj8npjPj3N0WGbjdlenK+1UPeXs/PfrjYtcXWY5Oc5yclxeW6eh6Qkj7H7qVAXZistUqgkj7JyEVmYvmiQTamfb3KsFAHuBoJemf1B28gDOdapq1N3Zt++YvC4A/KRSRYywY5q1XK+LthoBAE5WVvIH7fU3OVhZKfV6AEj3GpXuNYpZqdTrRz24P3p6enp6ejoAZGZmuri49Bukvr5+7dq1xcXF0dHRN27cSEhI8PPzmzFjRl5e3muvvfbWW28tXLgwLi5u1KhRA0UYHEVROp3Okv8Q09PT16xZc/78+aVLlx46dCg4ONjwUdKiRYvmz59fUlIya9Ys84Kr1WpUXH5bAAALsElEQVQ7O9OhfFhWVlYjRoxwcXHp7OwMCQlhKz2KovR6vbW1tYXpMUNcVFQUGBjIHE4MCQlhPpy3BCulY6Snpy9fvry8vHzZsmWHDx+eMmWKhQH1ej1N06I+/3Sah+lozBCzEpCV0rm7u1dUVLi4uCiVSrFYrNFo3NzcfjNDmqb7zoDEQoPOb++QPdiVWD7S3cFKwDxe4Oqc29pmw+dHjBjxcv31j1t+4fHAmsf/PNDfwUrw0iiPpKvV9lYCAY+X6OgQbmfnL7ZZLWv4tkPuIBC84e3VTdNmL57vVG1taqnt7l5YVfvcSPeFri4vXmv4S821borODQpcM9qTyXBmRfVLozxUFGW8LQBUqbv9xffb1trRnukNjYfutiv0+uzAgIuqro9abq8Z7WmIP04sru/pMbzWYjcXAJg4wm7HbdMbOV5WqVPc3YZe2OLiYqlUGh0dDQBjx46NjY0tLi6eMWPG7NmzZ8+eDQA0Tfv6+r7++utmj53Zqqurr1+//qc//UkoFE6fPj0uLq6wsNDNze3mzZuGfVKtVsvW7yEr2tra2tvbOZUeM8RRUVEqlYoZ4lOnTlneoFlhPMRTp06Ni4s7ceKE5Q36jyA4OPjjjz82LFZWVqalpZkZy/KPLFG/ZpRX9VAWncVRWlrq4ODAfKzf1tbm7e29Z8+eO3fu+Pr6XrlyRa/Xb968OTQ0lDJ6lYdiyVkcJSUlEonkwoULhBCZTObh4ZGfn3/8+HF7e/uysjJCyMmTJ21tbSsqKsyLT1j6PD0mJsZwFsfp06fZSo+tsziYIa6trZXL5cwQHzx40PKwrJTOMMRKpfLq1aseHh7ffvut5WFZPIuDEKLVaqOiorh2Fkdvb6+Xl9fnn3/e3t6+f/9+JycnuVz+m1v1exYHNuhHpfie8h83bhoW19RdL1f1M/aDNGhCSHZ2dlBQkL+//7hx4959993Ozs7u7u6cnBwfHx8HB4f4+PiGhgazM7SkQRNCdu7cGRgY6ODg4Ovru3XrVr1e39HRsXPnTn9/f0dHx7CwsPz8fLODE8t+VXp7e62tra2trXk8nlAotLa2fuedd3p6ethKj60GTQjJzs6WSqV+fn7MELMSk5UuQx4Msb29vY+Pz9atW1mJyVaD7jvEmzdvtjwsW6W7dOlSZGSkWCwOCgo6ffr0UDbpt0HzCCHs7NYjs2RmZq5bt26IT1apVFZWVsanjlnC8mPQJtGUSqWzszMr0YDVA6kA0NnZKRKJWHyzrByDZuj1epVK5eTkxEo0eASlY1ohK9FYPwbd1dXF2dIxx6CH+GZpms7Kynr55ZeNV+I8ogghxFHYoBFCiKOwQSOEEEdhg0YIIY7CBo0QQhyFDRohhDgKGzRCCHEUNmiEEOIobNAIIcRR2KARQoijsEEjhBBHPar5oNEQtbe3Z2ZmDvHJKpVKKBSyNaEEIYSmaYFAwEo01ufi0Ov1Vlas/XyyO6EEu6VjfS4O1kvH4jQmNE0TQtgqnU6nU6vVzO2fWMFu6R5qLg5CSD8/A6xM3YQej3/+858szqzIrhs3bsyaNWu4sxjQxo0bLZxd79GpqalJTk4e7iwGtG7duoKCguHOon+XL19esmTJcGcxoBdffPHHH3+0JAIe4kAIIY7CQxy/J9HR0Z6ensOdRf8kEslv3r5yGMXExPj6+g53Fv1zcnJKSkoa7iwGNGXKlDFjxgx3Fv1zdXWdOXPmcGcxoGnTpo0aNcqSCDgfNEIIcRQe4kAIIY7CQxy/D93d3ZmZmTweTyQSrVu3jq0Pwc1WVlaWm5vr4eEBAC+99JKdnZ1xehqNZliypWn6yJEj+fn5ubm50KdoJlk95iRNcuNaAZVK5ccffywSibRa7d///ncA4E7pTHKrqqriVOlaWlp27Nhha2vb09OzYcMGoVDIYukE77zzzqN+A8hyx44d8/f3T01NbWlpUalUw35MsK6uLjAwcOXKldOnT7ezszNJr7KycliyVSgU1tbWzPkk0KdoJlk95iRNcuNaAa9duzZu3LglS5a0tLRQFFVRUcGd0pnkplarOVU6uVw+ffr0OXPmXL9+3c7O7vLlyyyWDg9x/D40NTX5+/sDQEBAwM2bN4c7HVCr1SUlJe+//35WVpZOpzNJb7iydXJyCgkJMSwOntVjTtIkN64VcMKECZGRkXq9/vr164GBgZwqnUluXCudr6+vRqPZtGlTa2treHg4u6XDQxy/MzRN83i84c4CnnzyycmTJ7u6uubl5Z08edKw3iQ9jmQ7eFbDkiQHC6hUKnfs2JGSkmJ83QdHSmecGwdLN3r06C1btuzbt6+4uHiI+QwxPdyD/n3w8/Orr68HgPr6ej8/v+FOB27duqXRaABALBbTNG2SHkeyHTyr4U2SawVUqVSZmZnPP/98QEAAcKx0JrlxrXTffPNNZWUlADg5OanVanZLh6fZ/T709PR88sknNE3b2dmtXbuWzx/mv6w3b97csWOHRCLh8Xjp6ekAYJyeRqMZlmxlMtnhw4dramqCg4MTExMnTpw4SFaPOUmT3Dw9PTlVwH379l24cMHd3R0AZs6cGR4ezp3SmeQ2cuRITpWura3ts88+s7Gx0Wq1GzZs4PP5LJYOGzRCCHEUHuJACCGOwgaNEEIchQ0aIYQ4Chs0QghxFDZohBDiKGzQCD2E1NTUgIAAHx8foVAYEBAQEBCwY8eO4uLiv/zlL2bHpGk6OTn53LlzxivT0tJycnLMC/jqq69mZ2ebnQ/iDrySEKGHsHfvXgCQyWRJSUkNDQ2G9XFxcWbH/Oyzz4KCgqZMmcJCfgAAsHXr1okTJyYlJXl7e7MVEw0L3INGyFInT56cP39+cXHxnDlz1q9fn5iY+MILL7z//vsrVqyIiYlpb28HgM2bNyckJCQmJq5evVqn0xlvnpGR8eqrrwIARVHPPfdcbGzss88+e/v2bWbNiy++GBsb++STT65fvx4AYmNjDdcTJyUlFRQUtLW1zZ49Oz4+fsKECVlZWQAgFArT0tI+/vjjx1wHxDps0AixQyAQXL169f333y8sLDx48GBYWNiePXsmTpz4ww8/lJSUlJSUFBUVFRUViUSiPXv2GLaSyWSOjo7M5JkFBQW3bt06e/bs7t27q6urAaC9vV0qlZ49e/b8+fNnzpwpLy83HPqQy+V1dXWzZ88+fPhwSEjI6dOny8rKmGugAWDOnDknTpwYhiogVuEhDoRYExoaytw43NnZOSIiAgDc3NwUCkVJSUl1dXVCQgIAdHV1ubm5GTZpbm42HIioqKiYPHkyANjY2MTGxgKAi4tLR0fHokWLrK2tW1tb29vblyxZ8uabb6rV6ry8vNTUVIFAMHPmzIyMDKVSmZSUlJaWxoTy8fFpbm5+vO8esQ8bNEKsMZ5/3fCYEMLj8Z5++umMjIy+mzDfNTw2rKcoCgBycnLKy8u///57gUDAHKS2sbFZvHhxXl7eoUOHmAMa48aNq6mpKS4uzsvL27p168WLF7kwgyBiBR7iQOiRmzp1akFBQW9vLwBkZGRcvnzZ8C1vb++mpibmcUhISGlpKQCo1er//Oc/ANDW1ubv7y8QCK5cuVJbW6vVagHgr3/966effsrj8Zjp0A4fPlxeXj579uxPP/20qampq6sLAJqbm728vB77G0Uswz1ohB652NjY5cuXx8XFWVtbe3t7r1mzxvAtqVSqUCju3r3r5uaWnJx84MCBmJgYLy+vqVOn0jSdkpKyePHimTNnhoaGvvHGG5s2bZo0aZJUKqVpeuXKlUyEoKCgl156ycrKSqvVbty4USKRAEBhYeHs2bOH5c0iFuFsdggNs+3bt7e0tGzZsmWIz79x48a8efMqKyuFQmG/T9Dr9REREUePHvXx8WEvTTQM8BAHQsPs5Zdfrq6uNrlQZSDvvffe4sWLd+3aNVB3BoBNmzatW7cOu/P/B3APGiGEOAr3oBFCiKOwQSOEEEdhg0YIIY7CBo0QQhyFDRohhDjq/wH8Dkpuj2CCOwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# A legend is constructed\n",
    "\n",
    "legend_labels <- c(paste(gene, \"(+)\", sep=\" \"),\n",
    "                   paste(gene, \"(-)\", sep=\" \"))\n",
    "\n",
    "ggsurvplot(fit, data = all_clinical_df, pval = FALSE, legend = c(.75,.9), fun='pct', break.time.by=500, \n",
    "           risk.table = TRUE,\n",
    "          xlab = \"Time (days)\",\n",
    "          ggtheme = theme_bw(),\n",
    "          legend.labs = legend_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Call:\n",
       "coxph(formula = Surv(new_death, death_event) ~ event_rna, data = all_clinical_df)\n",
       "\n",
       "  n= 178, number of events= 58 \n",
       "\n",
       "                                   coef exp(coef) se(coef)      z Pr(>|z|)  \n",
       "event_rnalow_ENSG00000164078.11 -0.4484    0.6386   0.2689 -1.668   0.0954 .\n",
       "---\n",
       "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
       "\n",
       "                                exp(coef) exp(-coef) lower .95 upper .95\n",
       "event_rnalow_ENSG00000164078.11    0.6386      1.566     0.377     1.082\n",
       "\n",
       "Concordance= 0.556  (se = 0.036 )\n",
       "Likelihood ratio test= 2.83  on 1 df,   p=0.09\n",
       "Wald test            = 2.78  on 1 df,   p=0.1\n",
       "Score (logrank) test = 2.83  on 1 df,   p=0.09\n",
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# A cox regression can be performed to determine how significantly the expression patterns affect survival\n",
    "\n",
    "res.cox <- coxph(Surv(new_death, death_event) ~ event_rna, data=all_clinical_df)\n",
    "res.cox.extended <- summary(res.cox)\n",
    "res.cox.extended"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}