{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Regulome Explorer T test for numerical and binary data\n", "Check out more notebooks at our ['Regulome Explorer Repository'](https://github.com/isb-cgc/Community-Notebooks/tree/master/RegulomeExplorer)!\n", "\n", "In this notebook we describe the implementation of a student's t test to compute the significance of associations between a numerical feature (Gene expression, Somatic copy number, etc.) and Somatic mutation data which can be categorized into two groups according to the presence or absence of somatic mutations in a user defined gene. Detail of the test can be found the following link: https://en.wikipedia.org/wiki/Welch%27s_t-test\n", "\n", "To describe the implementation of the test using BigQuery, we will use Gene expresion data of a user defined gene and the precense or absence of somatic mutation of a user defined gene. This data is read from a BigQuery table in the pancancer-atlas dataset." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Authenticate with Google (IMPORTANT)\n", "The first step is to authorize access to BigQuery and the Google Cloud. For more information see ['Quick Start Guide to ISB-CGC'](https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowToGetStartedonISB-CGC.html) and alternative authentication methods can be found [here](https://googleapis.github.io/google-cloud-python/latest/core/auth.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Import Python libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from google.cloud import bigquery\n", "import numpy as np\n", "import pandas as pd\n", "from scipy import stats\n", "from scipy.stats import mstats\n", "import seaborn as sns\n", "import re_module.bq_functions as regulome" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## User defined Parameters\n", "The parameters for this experiment are the cancer type, the name of gene for which gene expression data will be obtained, and the clinical feature name. Categorical groups with number of samples smaller than 'MinSampleSize' will be ignored in the test. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cancer_type = 'LGG'\n", "gene_expre = 'DRG2'\n", "gene_mutation = 'TP53'\n", "MinSampleSize = 10\n", "\n", "bqclient = bigquery.Client()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data from BigQuery tables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step is to select all participants in the selected study with a least one mutation." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "barcode_set = \"\"\"barcodes AS (\n", " SELECT Tumor_SampleBarcode AS SampleBarcode \n", " FROM `pancancer-atlas.Filtered.MC3_MAF_V5_one_per_tumor_sample`\n", " WHERE Study = '{0}' \n", ")\n", "\"\"\".format(cancer_type)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Gene expression data from the BigQuery:** The following query string retrieves the gene expression data of the user specified gene ('gene_expre') from the 'Filtered.EBpp_AdjustPANCAN_IlluminaHiSeq_RNASeqV2_genExp_filtered' table available in pancancer-atlas dataset. The gene expression of a participant is computed as the average gene expression of the tumor samples of the participant. Moreover, we are considering only samples with a least somatic mutation. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "query_table1 = \"\"\"table1 AS (\n", "SELECT Symbol, data, ParticipantBarcode\n", "FROM ( \n", " SELECT \n", " Symbol AS symbol, AVG( LOG10( normalized_count + 1 )) AS data, ParticipantBarcode\n", " FROM `pancancer-atlas.Filtered.EBpp_AdjustPANCAN_IlluminaHiSeq_RNASeqV2_genExp_filtered` \n", " WHERE Study = '{0}' AND Symbol ='{1}' AND normalized_count IS NOT NULL\n", " AND SampleBarcode IN (SELECT * FROM barcodes)\n", " \n", " GROUP BY \n", " ParticipantBarcode, symbol\n", " )\n", ")\n", "\"\"\".format(cancer_type, gene_expre )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Somatic mutation data from the BigQuery:** The following string query will retrieve a table with patients with at least one Somatic mutation in the user defined gene ('mutation_name'). This information is extracted from the 'pancancer-atlas.Filtered.MC3_MAF_V5_one_per_tumor_sample' table, available in pancancer-atlas dataset. Notice that we only use samples in which FILTER = 'PASS'." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "query_table2 = \"\"\"table2 AS (\n", "SELECT \n", " ParticipantBarcode \n", "FROM\n", " (\n", " SELECT\n", " ParticipantBarcode AS ParticipantBarcode\n", " FROM `pancancer-atlas.Filtered.MC3_MAF_V5_one_per_tumor_sample`\n", " WHERE Study = '{0}' AND Hugo_Symbol = '{1}'\n", " AND FILTER = 'PASS' \n", " GROUP BY ParticipantBarcode\n", " ) \n", ")\n", "\"\"\".format(cancer_type, gene_mutation )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point we can take a look at the combined data (Gene expression and Somatic mutation) by using a simple LEFT JOIN command. Participants with and without somatic mutations are labeled as 'YES' and 'NO' respectively. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " in runQuery ... \n", " the results for this query were previously cached \n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
data1data2ParticipantBarcode
12.986545NOTCGA-VW-A7QS
22.398801YESTCGA-DB-5280
32.513401YESTCGA-HT-7611
43.063967NOTCGA-DH-5141
52.549122NOTCGA-HT-A616
62.769918NOTCGA-HT-7694
72.527377YESTCGA-CS-6290
82.832638NOTCGA-HW-7495
92.373629YESTCGA-WY-A85C
\n", "
" ], "text/plain": [ " data1 data2 ParticipantBarcode\n", "1 2.986545 NO TCGA-VW-A7QS\n", "2 2.398801 YES TCGA-DB-5280\n", "3 2.513401 YES TCGA-HT-7611\n", "4 3.063967 NO TCGA-DH-5141\n", "5 2.549122 NO TCGA-HT-A616\n", "6 2.769918 NO TCGA-HT-7694\n", "7 2.527377 YES TCGA-CS-6290\n", "8 2.832638 NO TCGA-HW-7495\n", "9 2.373629 YES TCGA-WY-A85C" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sql_data = 'WITH\\n'+ barcode_set+','+ query_table1+','+query_table2 \n", "\n", "sql = (sql_data + \n", "\"\"\"\n", "SELECT \n", " n1.data as data1, \n", " IF( n2.ParticipantBarcode is null, 'NO', 'YES') as data2, \n", " n1.ParticipantBarcode\n", "FROM\n", " table1 n1 \n", "LEFT JOIN table2 n2 \n", "ON n1.ParticipantBarcode = n2.ParticipantBarcode\n", "\"\"\")\n", "\n", "df_data = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n", "df_data[1:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To visualize the gene expression data in both groups with or without somatic mutation, we can use a 'violinplot' plot:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0lPed5/v3V1Wlfd+Q0IIWFiGBhEC2wTZgsN2G4PaSieM4jpfOYudM7nQynb7TfXrOTbqT02emt3Tfe3oyjjvp2Jk4idfuxnYwtrEdTAzYEmhFiF0s2lcktJbqd/+okiyEQAjrqSrp+b7OqWNV1e+p+mJQfer3PL9FjDEopZRSACGBLkAppVTw0FBQSik1QUNBKaXUBA0FpZRSEzQUlFJKTdBQUEopNUFDQSml1AQNBaWUUhM0FJRSSk1wBrqA2UpOTjY5OTmBLkMppeaVioqKDmNMykzt5l0o5OTkUF5eHugylFJqXhGRxutpp6ePlFJKTdBQUEopNUFDQSml1AQNBaWUUhMsCwURCReRj0WkSkTqROSvrtH2CyJiRKTMqnqUUkrNzMrRR8PAVmNMv4i4gH0isssYc2ByIxGJAf4YOGhhLUoppa6DZT0F49Xvu+vy3abb5u2HwN8CQ1bVopRS6vpYek1BRBwiUgm0Ae8YYw5Oeb4UyDLGvGFlHerqdDtWpdRkloaCMWbMGLMGyARuFpFV48+JSAjwj8B3Z3odEXlKRMpFpLy9vd26gm2mpqaGH/zgB+j/U6XUOL+MPjLG9AAfANsmPRwDrAI+EJEzwHpg53QXm40xzxpjyowxZSkpM87SVtfp0KFDjIyM0NHREehSlFJBwsrRRykiEu/7OQK4Czg6/rwxptcYk2yMyTHG5AAHgPuMMbqGhZ+NjY0FugSlVJCwsqeQDrwvItXAJ3ivKbwhIj8QkfssfF81S0NDeo1fKeVl2ZBUY0w1UDrN49+7Svs7rKpFXdvAwECgS1BKBQmd0Wxjw8PDAPT398/QUillFxoKNtZ38SIAF33/VUopDQWb8ng89PrCoKenJ8DVKKWChYaCTfX19U2MOuru6gpwNUqpYKGhYFNdviDIio+j/9IlRkZGAlyRUioYaCjYVGdnJwB5SYmX3VdK2ZuGgk11dnYSIkJOQjyAzmpWSgEaCrbV0dFBQmQEiZGRgPYUlFJeGgo21dnRQUJEOKFOB9FhYRoKSilAQ8GWjDF0dXcTHxEBQHxEuIaCUgrQULCl/v5+RkdHiY8IByA+PFyHpSqlAA0FW+ru7gYgLtwbCnER4fT19eF2uwNZllIqCGgo2ND4DOaJUAgPwwC9vb0BrEopFQw0FGxo/MM/NjzM+9+wsMseV0rZl4aCDfX29hLqdBLm9K6cHhOuoaCU8tJQsKHe3l5iwkIn7sdoT0Ep5WPldpzhIvKxiFSJSJ2I/NU0bf5ERI6ISLWI7BGRJVbVoz7lDYWwifsuh4Nwl0tDQSllaU9hGNhqjCkB1gDbRGT9lDaHgTJjTDHwCvC3FtajfHp7ei7rKYD3uoKGglLKslAwXuNberl8NzOlzfvGmPG9IA8AmVbVo7xGR0fpv3RpYuTRuNjwMHp8Q1WVUvZl6TUFEXGISCXQBrxjjDl4jeZfA3Zd5XWeEpFyESlvb2+3olTbGF8yO8E3m3lcfEQ4XV1deDyeQJSl1GU8Hg+vvvoqNTU1gS7FdiwNBWPMmDFmDd4ewM0ismq6diLyFaAM+LurvM6zxpgyY0xZSkqKdQXbQFtbGwCJkZeHQlJkJKNut+7CpoJCf38/hw4d4je/+U2gS7Edv4w+Msb0AB8A26Y+JyJ3Af8duM8YM+yPeuzs/PnzOEJCSI6OuuzxtJjoieeVCjTd9ClwrBx9lCIi8b6fI4C7gKNT2pQCP8EbCG1W1aI+dfr0adJionGGXP5XnxIdhcvh4PTp0wGqTKlPDQ0NBboE27Kyp5AOvC8i1cAneK8pvCEiPxCR+3xt/g6IBl4WkUoR2WlhPbbX1dXFhQsXWJqcdMVzjpAQ8hITOFJXp9cVVMBpKASO06oXNsZUA6XTPP69ST/fZdX7qysdPHgQAVameq/LvHvsJAB3Lc8HoDAtlYaaI9TW1lJcXByoMpWiv79/5kbKEjqj2SZ6e3s5cOAAhWmpxPmWzG7r76dt0i/fsuQkkqOi2PPuu7piqgqoixcvTvw8ODgYwErsR0PBBjweDy+//DJiDLflXH3SuIiwOT+Hjs5O3nnnHT9WqNTluifNmenW+TN+paGwwBljePvttzl9+jR3LcsnYcpQ1KmWJidRmpHOvn37qK6u9lOVSl2uo6ODkBAHoPuH+5uGwgJmjOGtt97iww8/pDQjndXpi67ruK1L88mMj+Oll16isrLS4iqVupwxhpbWVpIXLwERWltbA12SrWgoLFBDQ0O89NJL7Nu3j7UZi7l7+VJE5LqOdTpCeKhkFVnxcbz88svs2bNHRyQpv+nt7WXg0iVik1KJionnwoULgS7JViwbfaQC5/z587z4m9/Q3d3Nprwc1i/Juu5AGBfqcPCF4iJ2Nxznvffe49SpU3zxi18kLi7OoqqV8jp37hwAMQmp9Pd0cfbcOTweDyEh+h3WH/T/8gIyMDDAzp07eeaZZxgdHODLa0vYkJM960AY53I4uLewgM+tXM6Fc+f4p3/8R/bu3asjk5SlTp48idPpIiY+ifiUNIYGB2lpaQl0WbahPYUFYHR0lPLycvbs2cPQ0BBrM9K5PXcJ4S7XnLz+6vQ0suLjeO/4KXbv3s0nn3zCtm3bWLlypX57U3PKGMOxY8eJS05HQkKIT8kA4Pjx4yxevDjA1dmDhsI8NjIywscff8y+Dz+kr7+f7IR47lpdSMqUdY3mQnxEBJ8vLuJUZxd7TpziV7/6FampKdxxxxZWr16t4aDmREtLC729PSzP866dGRYRSUxCMkeOHGHz5s0Brs4eNBTmoZ6eHsrLyzl44AADg4MsSYhnR2kx2fFxN3yq6HrlJSWSk5DA0bZ2Pmo8x0svvcS7777LrbfeSmlpKeFT9mlQajZqampAhKT07InHkhfncLqunO7ubhISEgJYnT1oKMwTHo+HEydOcPDgQRoaGjDGkJ+UyIbCFWTExfq1lpAQoTAtlZWLUjjW3snBs+d444032L17N8XFxdxyyy1kZGT4tSY1/3k8HqqqqklISSc0/NP5NCmZuZyuK6e6ulp7C36goRDk2traqKqqorKykp6eHiJDQ7klO5OSxWnER1x7IprVRIQVqcmsSE2m5WIfhy80U1V5mIqKCtLT0yktLWX16tXExvo3tNT81NjYSE9PNwVlqy97PCIqltikRRw6dIhNmzZZ3hu2Ow2FIHTx4kVqamqoPHyYpuZmBFiSmMDGogJWpCTjCMLz92mxMWyPjWHL0jzqWlqpbW3jt7/9Lbt27SIvL4/S0lIKCwsJCwsLdKkqSFVUVOB0ukhenHvFc2lLlnHs0D7OnTtHdnb2NEeruaKhECS6u7s5cuQIdXV1NDY2ApAWE8PWZXmsTE0lOiw0wBVen3CXk3VZGazLyqDz0gB1rW0caWrilZMncTocLF22jKKiIgoKCoiMjAx0uSpIDA0NUVtbS0pWPg7nlR9LKRm5nKw+QHl5uYaCxTQUAsQYQ3t7O/X19dTV1nKhqQnwbnZze+4SClJTSIqa3x+aSVGRbMrLYWPuEi5cvEhDWwfHzpzh6NGjhIiQm5c3ERA6Kc7eqqurGR0dJW3J8mmfd7pCScnIo7q6mh07dmiP00KWhYKIhAN7gTDf+7xijPn+lDZhwC+AdUAn8LAx5oxVNQXa2NgYjY2NHD16lPojR+jyrf6YHhvDHfm5LE9JnnHBuvlIRMiMiyMzLo6tS/No6evnWHsHDc1N7Dx5kp07d7J48WJWrlxJQUEB6enpet7YZsorKoiKSyAmIfmqbdJyltPSeIyamhrKysr8WJ29WNlTGAa2GmP6RcQF7BORXcaYA5PafA3oNsYsFZEvAX8DPGxhTX43MDDA8ePHOXr0KMeOHWNoaAhHSAjZ8XGULV9KfnIisTYaxikipMfGkB4bw6a8HDouDXCio5MTnV3s2bOHPXv2EBsbS0FBAStWrCAvL4/Q0Plx6kzdmPb2di6cP0/e6puv+WUgNjGVyOg4KisrNRQsZOXOawYY38HF5buZKc3uB/7S9/MrwD+LiPiOnZeMMbS2ttLQ0MDRo0c5d+4cxhgiQ0NZlpTA0uQ8chISCHU6Al1qwIkIKdFRpERHsSEnm0sjI5zs6OJkZyeHKyr4+OOPcTqd5OXlsWLFClasWKHj1BegqqoqECE1M++a7USElKx8Ttcfoqenh/j4eD9VaC+WXlMQEQdQASwF/pcx5uCUJhnAOQBjjFtEeoEkoMPKuuba6OgoJ0+epKGhgYajR+n17Rq1KCaaDUuyyE9KJD02Rk+JzCAqNJTixWkUL07D7fFwrruHk51dnDx3jmPHjvH666+TmpJCge80U1ZWls6kXgBq6+qIS1pEWMTMM/FTM3NprD9EfX09GzZs8EN19mNpKBhjxoA1IhIP/JuIrDLG1E5qMt2n5BW9BBF5CngKCJqRB729vRO9gVMnTzLqdhPqcJCTGM+GgmXkJyUSrRfDbpgzJITcpERykxK50xi6BgY52dnFqc4u9n34IXv37iUyIoJly5ezcuVKli1bprOp56Guri7a29rIL77lutpHxsQTGROvoWAhv4w+Msb0iMgHwDZgciicB7KA8yLiBOKArmmOfxZ4FqCsrCxgp5ba29snho2Or/EeFxFOcVoq+UmJZCXE45wn31zfPXaStj7v2b1fHaoiNTqau5bnB7iq6YkISVGRJEVFcnN2JsNuN6e7ujnR0cmxI0eoqqoiJCSE3NxcioqKKCwsJCYmJtBlq+tw8uRJABJSM6/7mITUxTQ2HsPtduOcZviq+mysHH2UAoz6AiECuAvvheTJdgJPAPuBLwDvBdP1BGMMLS0t1NbWcqSujrb2doCJi6RLk5NIjoqcl6eF2vr7GR4bA+BcT2+Aq5mdMKeTgtQUClJT8BhD08WLnGjv5Nik0UzZ2dkUFRVRVFSk1yGC2NmzZwkNCycy5vqHJMclp3Ph5BGam5vJysqysDp7sjJm04HnfdcVQoCXjDFviMgPgHJjzE7gZ8D/EZETeHsIX7KwnuvW19dHVVUVhw5V0NrahgBZ8XHctSyfZSlJthotFOxCJg133ZyfS8elAY61d3CsvZNdu3axa9cucnNzWbt2LUVFRTq+Pcg0NTURFZc0qy9WMfFJE8dqKMw9K0cfVQOl0zz+vUk/DwEPWVXDbHg8HhoaGvjkk084fuwYHmNIj43hD5YvZUVqMpE6LDLoTR7NdFvuEroHBqlvbaOmpYVXX32V13fupGjVKm655Rb9MAkCxhg6OjpYlLNiVseFRUbjcDjp6JhX41HmDdufkPN4PNTW1vLB++/T2tZGTFgYN2dnsipt0byfUWx3CZER3Jq7hA052VzovUhtSytHamo4fPgw+fn5bNmyhdzcK9fZUf4xMDCA2+0mPHJ2139EhPDIaHp6eiyqzN5sHQonT57k9Z07ae/oICkqknsLV7AyNZWQkPl3jUBdnYiQGR9HZnwcW5fmc7ipiU/OnuOnP/0pOTk5PPjggyQnX30mrbLGwMAAAKFhsz8d6wwLnzhezS1bhoIxht///ve89dZbJEREcP+qlaxISZ6XF4zV7IQ6HdySncXajMVUNbXw0Zmz/PjHP+bhhx9mxYrZncZQn83w8DAAIc4rt409UeVd+GBpyfppj3U4XQwPj1hXnI3Nj/GTc+ydd95h165dLEtO4omyUgpSUzQQbMblcFCWlcETZaXEh7r4xS9+QW1t7cwHqjkzPtBwut+9/t5O+ns7r3qsiGCMx7La7Mx2oeB2uzl44ADLU5J5YNVKXW7C5uIiwnl0bQmJkZHs378/0OXYykQY3MAodGOMfpGziO1C4cKFCwwND7MseXbD4BaaYbeb0NBQNmzYQGhoKMNud6BLChiXw0FuYjyNZ84wOjoa6HJsw+XynjYaG5v9vz3PmFsXSrSI7UIhLS2N6OhoDl1oIojmyfndsNvNunXruPfee1m3bp2tQ2FwdJS61naW5OToDFk/Gp8zMjY6+2sDY6OjOufEIrYLhbCwMLZv307zxT5erKzh4tBwoEsKiDCnk4qKCt544w0qKioIs+mHYfPFPn5ZUcXI2Bh/+Id/aOveo7+N77w3OjL730H36DARAd6jfKGyXSgAlJSU8MADD9Dcf4mff1JBdVMLHo+9eg1hTicjIyPs37+fkZER24XCsNvNvlNn+GVFJW6nkyeffJK0tLRAl2UrLpcLp9N5Q6EwMjxEVNTMq6qq2bPXJ4GPiHDTTTeRl5fHyy+/zK6jx/io8RzrszNYlZaG02HLrLSFwdFRKs5doOJCM0OjoxQXF3Pffffpt84AEBGioqIYHR6c1XFjY27G3KMaChaxZSiMS0pK4qmnnqKhoYEPPviA3Q0n+P2Zc6xZnMbq9DRiw/Wc5ULR3n+JqqZmqptbGR0bo7CwkM2bN5OZef2rc6q5Fx0dzcDw0KyOGQ+R6OhoK0qyPVuHAkBISMjE3sCnT5/md7/7HftOnOD3pxvJS0qkeHEa+UmJOObJktjqU8NuN/Wt7VQ3t9B8sQ+Hw8GqVavYvHkzixYtCnR5CoiKiqK37erzEaYz6pv0Nn5NQs0t24fCOBEhLy+PvLw8urq6qKiooKKign+rOUJkaCgrU5MpSltEWky0XowMYh6P4Ux3N0da2jjW0cno2BipKSl87nMbKS0t1Q+SIBMZGcnYaPOsjhm/BqF/l9bQUJhGYmIid999N1u3buX48eNUVFRQ2dBAxfkmEiIjKFyUStGiVBIi9Tx0MDDG0NzXx5GWNurbOhgYGSE8PJyS0lLWrVtHVlaWBnmQCgsLw+2e3dyQMbd3CKvutGcNDYVrcDgcFBQUUFBQwODgILW1tVRVVfH706f5/elG0mNjJjZ70esP/tfef4n6Vm8Q9AwOTvx9lZSUsGLFCp1zMA+4XC7GZjlHZrz9+OQ3Nbf0t+Y6RUREcNNNN3HTTTfR09NDTU0N1dXVvH/iFO+fOEVmfByFqSm694LFugcGqW9rp761nY5LlyZO+20pLqaoqEhHEc0z3h7cjQ0H196fNazcjjML+AWQBniAZ40x/++UNnHAL4FsXy1/b4z5uVU1zZX4+Hg2btzIxo0b6ejooLq6muqqKt4+doJ3jp8kJyGelYtSWJ6SHLTj/1Ojoyf2aE6NiSY1iEdy9A0PU9/aztG2dpov9gGwJDubDVu3smrVKh2FMo95PB6E2X24j4eBx6ML4lnByk8sN/BdY8whEYkBKkTkHWPMkUltvgUcMcb8oW9P5wYRecEYM2/WxE1OTmbr1q1s2bKF1tbWiYD4bf0xdjecIC8xgcK0VPKTEnE5gmfxvbuW59PW7w2FL68tCXA1VxocHaWhrZ0jre0Te0gvXryYbbfexurVq4mPjw9whWouDA0N4Zxlz9rhCp04Vs09K7fjbAaafT/3iUg9kAFMDgUDxIg3+qPx7tM8LxfhERHS0tJIS0vj7rvv5vz581RXV1NTXc3x2npCHQ6WpSRRuCiVnIQE3chnGiPuMU50dHKktY3TXd14jCElOZk777yT4uJi3QhnAbp06RKu0NldjxvflKff96VGzS2/nNsQkRy8+zUfnPLUPwM7gSYgBnjYTLNIuog8BTwFkJ2dbWWpc0JEyMrKIisri+3bt3P69Gmqqqqoq62lrqWNyNBQClKTKVqUSnpsjK3PjXo8htPd3RxpaeV4RxejY2PExcZy2+23U1JSQlpamq3//yx0HR0dhEfFzuqY8Cjv9p2dnbOb36Cuj+WhICLRwKvAd4wxF6c8fQ9QCWwF8oF3ROTDqe2MMc8CzwKUlZXNq0WKQkJCyM/PJz8/n/vuu49jx45RVVVFdX09h2w6xNUYQ0tfP3UtrRNDSCPCw1mzdi1r1qwhOzubEJ0suOANDw/T2dlJ1orZzSp3hYYTGhZBc/Ps5jeo62NpKIiIC28gvGCMeW2aJn8E/E/jXcP6hIicBgqAj62sK1CcTieFhYUUFhYyNDREbW0tlZWVE0NcF8fFsjotlYLUVMJdwXmB+rO4ODREXUsbtS1tdA0M4HQ4WFFQwJo1a1i+fLkOIbWZs2fPYowhLml2s8tFhNikRZw6dVo327GAlaOPBPgZUG+M+dFVmp0F7gQ+FJFFwArglFU1BZPw8HDKysooKyujp6eH6upqDh86xO6GE7x7/BRLkxNZnbaI3MTEeX39YWRsjGNtHdS2tNLY3QNATk4Om9asYdWqVTqE1Mbq6upwOJ3EJc1+ddrERRkcO3yG5uZmFi9ebEF19mXlV7PbgMeAGhGp9D32F3iHn2KMeQb4IfCciNQAAvyZMabDwpqCUnx8PJs2bWLjxo1cuHCByspKqioraWirIzosjNVpqaxOT5s3p5fGTw9VNTVT39bBiNtNQkICW7dupbS0lMTExECXqAJsZGSE2tpaEtOycNxADzF5cQ7HK/dz+PBhDYU5ZuXoo31w7QHIxpgm4A+sqmG+EREyMzPJzMxk27ZtNDQ0UFFRwYFjx9jfeI7shHhK0tNYkZoclAv0Dbvd1DS3Ut3cQnv/JVxOJ6tWr6asrIwlS5ZoN19NOHz4MIODgyzPW3lDx7vCwknJzKW8vJw777xTl7yYQ3oSN0g5nU6KioooKiqit7eXQ4cOUVFezutHjvLeiVBKFqexJiOdmCDYkrC9/xKHLjRR19LG6NgYGYsXc/+dd1FcXKy/rOoKIyMjfPDB74hJSLmhU0fjMpeuou3cST788EPuvvvuOazQ3jQU5oG4uDi2bNnC5s2bOXnyJPv37+ejhgb2N55jeUoS67OzSIuN8WtNxhhOdXbx8dnznO3pxel0UlxSwvr168nIyPBrLWp++fDDD7l4sZeSTZ/7TL3HmIRkUrPy2bdvH2VlZSQkJMxhlfaloTCPhISEsGzZMpYtW0ZXVxcHDx6k/JNPaCg/TG5iAhuWZJGVYO1MX48xNLS1c6DxPG39/cTFxbFt2zbWrVunSxmrGTU1NfHBBx+QkplHfHL6Z3693KIyuprP8sorr/C1r31NhzLPAQ2FeSoxMZHt27ezZcsWDh48yO/37eNXh6vJio9j69K8Oe85GGM40dHJByfP0DUwQEpyMl/Yto3i4mIcQbR8hwpew8PDvPTSS7hCw1lWsmFOXjM8Mpr84vU0HPqQvXv3cscdd8zJ69qZhsI8Fx4ezubNm9mwYQMVFRW8/957PF9+mNXpi9iUl0P0HFxzaOvv573jp2js7iElOZlH7r+fwsJC/VamrpvH4+GVV16hvaOD4tu24Qqbu2tNi5Yso7vtAu+88w5paWkUFBTM2WvbkYbCAhEaGsqGDRsoLS3lgw8+4KOPPqKhvZN7li+lMC31hl7T4zF8dKaRj86cJTwinHvvvZebb75ZewZqVowx7N69myNHjpBffAsJqXM7hFREWL52IwP9vbz44ot8/etf1+tan4F+1VtgwsPD2bZtG9/+9rdJW7yY148c5c0jRxlxj83qdS4ODfHrw1X8/sxZStas4U/+5Lts2LBBA0HN2t69e9m3bx+L8wrJyC+y5D0cTierNtxNiCuM5557jvb2dkvexw40FBaopKQkvv71r7NlyxbqWtt5saqG4evc4aprYIBfVlTRNjjEQw89xEMPPaQXkdUN2bt3L2+//TapWfksLVlv6VyVsIgoVt92D26P4ac//SltbW2WvddCpqGwgDkcDu666y4eeeQRWvr6ebFy5mDoGhjk14dr8DgcPP3006xZs8ZP1aqFxBjDe++9x+7du0nNzKNg3Sa/TF6MjI6j+PbtjLjH+OlPf0pLS4vl77nQXDMURCRWRP6HiPwfEfnylOd+bG1paq4UFRXx5S9/mZa+fvYcPznxeGr05TuujXk87KyrxxMSwte+/nXS0m58YpGyL4/Hw5tvvsmePXtYlL2Ugps2I34clBAVm0DJxs/h9sC//Mu/cObMGb+990Iw09/Uz/EuVfEq8CUReVVExoezrLe0MjWnVq5cyaZNm6hpbuVUZxfg3X3truX5E20+Pnue1r5+7n/gARYtmt3KlUoBuN1uXn75Zfbv30/G0iJWrNuEiP9PSETGxLNm8w5CnGH8/Oc/p76+3u81zFcz/W3lG2P+3Bjz78aY+4BDwHsikuSH2tQc27p1K3FxcXxy7sIVz415PJSfb2LFihWsWrUqANWp+W5wcJDnnnuO6upqcovKyF99S0DXuwqPjKFk8w4iYhJ44YUXOHhw6h5fajozhUKYTIp5Y8xf493sZi+gwTDPOJ1O1q5dy5mubvqHhy977kxXDwMjI9x0000Bqk7NZz09PTz77LOcaWykoGwz2StKgmIBxNCwCEo2bidhUSY7d+5k9+7deDxXbO6oJpkpFF7HuyvaBGPM88B3gRGrilLWyc/3ni5q77902eNtvv1ux59X6no1NTXxzDPP0NXdw+pb72FR9tJAl3QZh9PFqvV3kZ5bwN69e3nppZdwX+dIPDu65uQ1Y8x/u8rjbwHLLKlIWSo+3rs2Uu/Q5T2Fi0PDREVGEhoaGoiy1Dx14sQJXnjhBcThYs2mHUTFBedeGRISwrI1txIeGU1NTTl9fX185Stf0U2epjHT6KNMEbl90v0/EZHv+W7X/DogIlki8r6I1ItInYh8+yrt7hCRSl+b393YH0PNVsiUrr2IdxihUtfr8OHDPP/884RGRFN6xx8GbSCMExGyV5RQcNMdNJ49y0+efZbe3t5AlxV0Zjp99HfA5GU3nwYuAQb4qxmOdQPfNcasxDtS6VsiUji5gYjEAz8G7jPGFAEPzaJ2dQMGBgYACJuy21W408nQ8DBjY7Ob+azsad++fbzyyivEJi2iZNMOwiKi5uy1T1QdoL+nk/6eTir3vsmJqgNz9toAi7LyWX3bPXR39/DMMz/R2c9TzBQKK4wxb0y6P2CM+QdjzA/xbat5NcaYZmPMId/PfUA9MHVBki8Drxljzvra6RREi3V2dgJDm58vAAAa1ElEQVSQEHn5gmTxERF4PB56enoCUZaaJ4wxvP322+zatYvkjBxW33oPTtfcnnLs7+1kzD3KmHuU3o4W+ns75/T1ARJSFlOy8XMMjYzwk5/8hAsXrhyRZ1czhcLUpQzvnPTzdY8+EpEcoBSYOiZsOZAgIh+ISIWIPH69r6luTFtbGwIkTDmXmhQVOfG8UtMxxvDmm2/yu9/9jvScFRTevIWQebwWVnR8Ems23YsRBz/72c84e/ZsoEsKCjOFQp+ILB+/Y4zpAhCRAqD/et5ARKLxTn77jjHm4pSnncA6YAdwD/D/TH6/Sa/xlIiUi0i5dvU+m5aWFhKjInFN+WVO9oWCLgugpmOM4fXXX5+YlLas9LaATEqbaxHRsZRs2kGIK4x//dd/pbGxMdAlBdxMf6vfB94QkSdEZLXv9iSw0/fcNYmIC28gvGCMeW2aJueBt4wxl4wxHXjnP5RMbWSMedYYU2aMKUtJSZnpbdU1dHR0kDjNiIswp5OY8PCJ00tKjRvvIRw8eJDMZasDPiltroVHRlOycQeusEiee+452/cYrhkKvqGnn8d72ug5320L8HljzK5rHSvefzU/A+qNMT+6SrP/ADaKiFNEIoFb8F57UBbp7ekhNnz6DU5iw0Lp7u72c0Uq2O3Zs2eih5C36qYFFQjjwiIiKd64HYcrnOeff97WPeYZ+3/GmFpjzOPGmHW+2xPGmFoRWTLDobcBjwFbfUNOK0XkcyLyTRH5pu+164G3gGrgY+Cnxpjaz/hnUldhjGFkdJQw5/TngcOcTkZGdE6i+tT+/ft5//33SVuyfMH1EKYKi4ii+PZtGAnh5z9/zrZfkGbceU1ENuAdNbTXGNMmIsXAnwMbgayrHWeM2Yd3Mb1rMsb8Hd6hr8pi4/MQrvaXIuhcBfWp+vp63njjDZLSs1leetuCDoRx4VExrL71Hir3vsnzzz/P008/bbsJbjNNXvs74F+B/wS8KSLfB97BO4pIZzTPMyEhIURGRNB/ld5A/8gI0ZOW0lb21dzczIsvvkhMQgorb9ri16WvAy0qLpHCW+6ko6OTX//617ZbK2mmnsIOoNQYMyQiCUATUGyMOW59acoKCYmJdF+6cuCYMYaewSGWJCQEoCoVTAYGBvjlL18gxOmiaP1dOJz228o9IXUxy9bcyrHD+3jnnXe45557Al2S38wU/4PGmCEAY0w30KCBML9lZ2fT3NfP2JRvPx2XBhh2u8nOvuacRLXAGWN45ZVX6L3YS+EtdxIWYd9tWNNzV0wsomen/Rhm3E9BRHb6bq8DOZPu7/RHgWpu5ebmMjo2xoXey6eMnO7yXlTLyckJQFUqWBw4cICGhgbyVt1EbGJqoMsJuKXF64mOT+K1117j4sWp06wWppn6hfdPuf/3VhWi/GPp0qU4HA6Od3SSnfDpslYnOjpJW7SIBD19ZFtdXV289dZbJC7KJCO/KNDlBIUQh4OVN91BxXv/wc6dO/nKV74S6JIsN9PS2ROrlopIiu8xnVI8j4WFhZGfn8/xc+fYujQPEeHSyAjne3q5Y8vaQJenAsQYw3/s3IlBWF56uy1GGl2vyJh4lqwspb72E+rq6igqWtiBOdPoIxGR74tIB3AUOCYi7SLyPf+Up6xQVFRE7+Agbb6Ndk50dGKAwsLCax+oFqyTJ09y4vhxclauJSxy7lY8XSgyl64iKjaBt3bvXvArCc90TeE7wO3ATcaYJGNMAt5Zx7eJyH+1vDpliRUrVgBwqrNr4r+xsbGkp6cHsiwVQO+8+y7hkdEszlsZ6FKCUkhICDmF6+jq7KSysjLQ5VhqplB4HHjEGHN6/AFjzCngK77n1DwUExNDeno6Z7q6McbQ2N3L8uXL9ZSBTTU1NXH+3Dkylq6a16ueWi0pPZuo2AQOHvw40KVYaqZQcPkWqruM77qCy5qSlD9kZ2fT0tc/MRR1yZKZVi1RC1VVVRUhISFBt7dysBER0pYs58KF83R0XPGxuGDMFArXWghHF8mZxxYvXszI2BjHOzon7it7On36NDGJqbhCwwJdCgDu0RFCQ0PZsGEDoaGhuEeD56MmMS0TgDNnzgS2EAvNNCS1RESmG5wrXLkBj5pHEhO9++me7fbutKZDUe2rtbWVtJyCQJcxwT06wrp167j33nsBOFRVE+CKPhURHYfD6aS1tTXQpVhmpiGpeoJxgYqJiQG8M5lDXS7CwoLjW6Lyr7GxMdxuN87Qud1S87NwukKpqKgAoKKiAldE8KzHJSI4XWEMDQ0FuhTL2GeVK3UZl8t7SejSiLerruwrJCSEMXfwDLN0ukIZGRlh//79jIyMzPke0J+VZ8yNcwGvB6WhYFM60kgBOBwOEhISGOjrCXQp88LI0CCjI8MkJV33FvXzjoaCTY1PwAl1OBb8ZBx1bbm5ufS0NzHmdge6lKDX2ezdqnMhrxFmWSiISJaIvC8i9SJSJyLfvkbbm0RkTES+YFU96nLj50RjwsIYHh7WzXVsrKSkhDH3KG3nTwW6lKBmjKH5TAOJSUlkZGQEuhzLWNlTcAPfNcasBNYD3xKRK9ZREBEH8DfAbgtrUVOMh0JcRDgeYxgeHg5wRSpQcnNzycjI4OzRQ3jGtLdwNZ1NjfR1t7Np48YFffrVslAwxjQbYw75fu4D6vFu6znVfwFeBdqsqkVd6dIl77pHSZHe9fIHBgYCWY4KIBFh27ZtDA1c4vSRQ4EuJyiNjgxzovoAKSkprF27sBeO9Ms1BRHJAUrxbuM5+fEM4EHgmRmOf0pEykWkvL1dF2mdC+OhkBwVedl9ZU95eXncfPPNnD9eQ1fr+UCXE1SMMRw//HtGhwf5whe+gGOBLwVieSiISDTensB3jDFTJ8L9E/BnxphrXuk0xjxrjCkzxpSlpKRYVaqtjIdAYlTEZfeVfW3fvp2U1FTqP35fRyNN0ni0kvYLp7n77rvJzMwMdDmWszQURMSFNxBeMMa8Nk2TMuA3InIG+ALwYxF5wMqalNfAwADhLhdRvjHgevpIhYaG8sTjjxMW6qL2o7cZHtQvCi2Nx2msP8SaNWvYuHFjoMvxCytHHwnwM6DeGPOj6doYY3KNMTnGmBzgFeA/G2P+3aqa1KcGBwcJczqJcDkn7iuVkJDAY489xtjoMNX73mJk2L7/LtrPn6bh0Ifk5eXx4IMPLuiLy5NZ2VO4DXgM2Coilb7b50TkmyLyTQvfV12HkZERQh0huHznR0dGgmfRMRVYWVlZPP7444wM9lO9bxcjQ/YLhvYLp6n/5AOys7N57LHHFvQM5qks+5MaY/bhXTjvets/aVUt6kputxtnSAghIojvvlLjcnNzefzxx/nFL35B1Ye/pfj2bYRF2GNHttazJzhasZfsrCyeePxx2y0DozOabU5EwCbdYjU7+fn5PPnkk7iHB6ja+yaDl6ZbMHnuRccl4XC6cDhdxCWnER3nvyUlmk7Vc7T8d+Tm5PDkk08SHm6/xaA1FGxKRDDGO9zOGGOb86VqdnJzc/nqV78KHjdVv3uTSxe7LX/PpSXriY5PIjo+iTWbdrC0ZL3l72mM4WxDFccrP2LFihU88cQTtl05WEPBpkJDQxn1jDHq8QDY9hdAzSwrK4tvfOMbuJwOqva+ycWuhTXP1BjDqdpPOF1XTnFxMY8++ujEKsJ2pKFgU+Hh4Qy5xxj2XUvQUFDXsmjRIp5++imio6Ko/nAXXa0XAl3SnDAeD8cO7eP88RrWr1/PQw89tOAnp81EQ8GmoqKiGBgZ4ZJv1FFUlD0uIqobl5iYyNNPP0VychK1+9+m/cLpQJf0mXjGxjjy8fu0NB5jy5Yt3HvvvYSE6Eei/h+wqejoaIwxdPR7J62N78Sm1LXExMTwjW98g8yMDN8H6vFAl3RDxtxuave/Q0fTGXbs2MFdd92l19V8NBRsKjY2FoCWvj5AQ0Fdv4iICL761a+Sn5dHQ8Vemk7VB7qkWXGPjlDz0W562pt48MEHufXWWwNdUlDRULCp8RBo6eu/7L5S1yM0NJTHHnuMFQUFHK/8iPMn6gJd0nVxj45Q8/vdXOxs5Ytf/CJlZWWBLinoaCjY1HgItPb1ExEebqsZm2puuFwuvvzIIxQWFnKy+gDnT9QGuqRr8gbCW/T3dPDII49QXFwc6JKCkoaCTY1fWHZ7PETqRWZ1g5xOJ1/60pcoKiriZPVBLpw8EuiSpjV+yqi/p5NHHnmEoqKiQJcUtDQUbCo0NJQQ34W1CBvO2lRzx+Fw8PDDD1NQUMCJqv20nDkW6JIuMzbmpu7Au/R1tfPwww9TWHjFBpBqEg0FmxKRiVNGLput7aLmnsPh4JFHHmHpsmU0HN5H+4UzgS4JAI/HQ/3B9+lpb+ahhx5i1apVgS4p6Gko2Nj4mGwdm63mgtPp5NEvf5mszEyOfuL9IA6k8R3TOlvOct9991FSUhLQeuYL/TSwMY/xLnHh8S11odRnFRoayuOPP05iUhJ1B971y1pJV9N49PDExLRbbrklYHXMNxoKNmWMwT3qXeJidHQ0wNWohSQyMpInn3iCsFAXdfvfYXR4yO81tJ47SWP9YUpLS7nzzjv9/v7zmZU7r2WJyPsiUi8idSLy7WnaPCoi1b7bRyKi/Ts/cbvdeIwBYHjI/7+0amEb38FtZGiAuoN7/NobvdjVzrFDH7JkyRIeeOABnak8S1b2FNzAd40xK4H1wLdEZOpl/9PAZmNMMfBD4FkL61GTTN6TWbfiVFbIysri85//PL0dLZyq/dgv7zk6PET9x3uIiY7h0Ucf1fk3N8CyUDDGNBtjDvl+7gPqgYwpbT4yxoyfdDwAZFpVj7rceBDER4QzODiI8fUalJpLa9asYf369Vw4UWf5AnrGGOo/+YDR4SEeffTLusjjDfLLNQURyQFKgYPXaPY1YJc/6lGf9hQSIiIY83h0j2Zlme3bt5ORkcGxQ/sYGui37H3On6ilu+0CO3bsICMjY+YD1LQsDwURiQZeBb5jjJl2Pz8R2YI3FP7sKs8/JSLlIlLe3t5uXbE2MrmnMPm+UnPN6XTy8MMPEyLQULHXkl5pf28Xp+vKKSws5Oabb57z17cTS0NBRFx4A+EFY8xrV2lTDPwUuN8Y0zldG2PMs8aYMmNMWUpKinUF28jw8DAAMb7NdcbvK2WFpKQkduzYQU9785yvqurxeGio2EtkRAQPPvigXlj+jKwcfSTAz4B6Y8yPrtImG3gNeMwYE1xz4xe48WGoEb5tB3VYqrLaunXrWLpsGafrPpnT00jnj9fQ39PJ/fffT2Rk5Jy9rl1Z2VO4DXgM2Coilb7b50TkmyLyTV+b7wFJwI99z5dbWI+aZHyIoMsRctl9pawiIjxw//0IcKJq/5y85uClPhqPVlJYWKiL3M0Ry8ZrGWP2Adfsxxljvg583aoa1NWNd7HHT+9ql1v5Q0JCAnfeeSe7d++mq/U8iYs+24DDUzUf4wgRduzYMUcVKp3RbFPj6x2Njo1ddl8pq916660kJCRyquYg5jP0UHvam+loOsPmzZuJj4+fwwrtTT8JbMrlu5Yw4LuWMH5fKas5nU62b9/GpYs9tJy9sT2ejTGcqv2E2NhYbr/99jmu0N40FGwq1Ldc9iXf/IRQXT5b+VFhYSGZmZk01h/G4+utzkZnUyN93e3ceeed+oVmjmko2FS4b2Od/uGRy+4r5Q8iwt13383w4CWazzTM6lhjDI1HD5OYlERpaalFFdqXhoJNhfnmJ1z0zU/QnoLyt/z8fLKzszl/vBqP5/p7C53NZ+nv7WLrli04HA4LK7QnDQWbGg+F/uFhQkNdeqFZ+Z2IcMcddzA0cIm2cyev6xhjDOeOVRMfH09xcbHFFdqTfhLY1HjPYGBklFCX9hJUYCxfvpxFixZx/kTtdS1/cbGrjYtdbdx+++3aS7CIhoJNjV+cM6DLC6uAERFuvfVWLvV2X9f2nRdO1BEeHs7atWv9UJ09aSjY1ORvWRoKKpBKSkqIiIiYcU2k4cEBOpoaWbdu3cTpTzX3NBRsanIohGg3XAWQy+Vi3bp1dDY3MjI0cNV2rWePY4xHV0G1mIaCTU1e1iJEl7hQAbZu3TqMMbRe5YKzMYbWxuPk5OSQnJzs5+rsRUPBpiaPNtKRRyrQUlNTycjIuOoopL7uDgb6e3Vegh/op4FNichEb0FPH6lgUFxcTH9PJ4P9V+7F1X7hNCEhIboSqh9oKNiYw9dD0KF9KhiMf+B3Np+94rmulrPk5uYSERHh77JsR0PBxsZHHWkoqGCQkJBAckoKXa3niI5LIjouCfDumTDQ10tBQUGAK7QHK3deyxKR90WkXkTqROTb07QREfn/ROSEiFSLiA4+9qPxUNAFxVSwWLZ0Kb2dbeStvomlJesBJuYvLF26NJCl2YaVPQU38F1jzEpgPfAtESmc0mY7sMx3ewr43xbWo6YI9Y311nWPVLDIycnBM+amv+fT7dp7O1qIiIxE92f3D8tCwRjTbIw55Pu5D6gHMqY0ux/4hfE6AMSLSLpVNanLhfnCQCcCqWCRmendia2vu2PisUu9HWRlZurugH7il2sKIpIDlAIHpzyVAZybdP88VwaHskiYb7lsDQUVLOLi4oiIiOBSbxcAHs8Yly72kJ6u3xX9xfJQEJFo4FXgO8aYqWPNpov+K1bFEpGnRKRcRMrb29utKNOWxk8b6V4KKliICKmpqQz09QAw2H8RYwypqakBrsw+LA0FEXHhDYQXjDGvTdPkPJA16X4m0DS1kTHmWWNMmTGmTM8rzp3x7rj2FFQwSUpKYuhSH+AdeQSQmJgYyJJsxcrRRwL8DKg3xvzoKs12Ao/7RiGtB3qNMTMvlajmxPhSxRoKKpjEx8czPDSAxzPG8ED/xGPKP6xcHvM24DGgRkQqfY/9BZANYIx5Bvgt8DngBDAA/JGF9air0FBQwSQmJgaAkaFBRoYGEBGio6MDXJV9WBYKxph9TH/NYHIbA3zLqhrU9dGls1UwGQ+A0eFBRoeHiIyK0vW5/Ej/Tysd6qeCSlRUFACjI8OMjgwTFRkZ4IrsRUPBxsbHhEfqL50KIuPrG7l9oaDrHfmXnjewsU2bNpGfn09WVtbMjZXyk/Eh0u7REcbcI0RExAW4InvRULAxl8tFTk5OoMtQ6jKXhcLoiA6E8DM9faSUCioulwsRwT06int0RCdX+pmGglIqqIgIYWFhjLlHcI+Oaij4mYaCUirohIWHMzI0iDEeDQU/01BQSgWdiPBwhnyzmTUU/EtDQSkVdCIiIiaWuNBQ8C8NBaVU0ImIiGB0ZGjiZ+U/GgpKqaAzuXegPQX/0lBQSgWdyb0DDQX/0lBQSgWdyUGgp4/8S0NBKRV0tKcQOBoKSqmgMzkIdGl3/9JQUEoFHe0dBI6V23H+q4i0iUjtVZ6PE5HXRaRKROpERHddU0oBuhtgIFnZU3gO2HaN578FHDHGlAB3AP8gIqEW1qOUmidCQ/WjIFAsCwVjzF6g61pNgBjxbvsV7WvrtqoepdT8oaEQOIG8gvPPwE6gCYgBHjbGeAJYj1IqSLhcrkCXYFuBvNB8D1AJLAbWAP8sIrHTNRSRp0SkXETK29vb/VmjUioANBQCJ5Ch8EfAa8brBHAaKJiuoTHmWWNMmTGmLCUlxa9FKqX8z+FwBLoE2wpkKJwF7gQQkUXACuBUAOtRSgWJ8Z5CcXFxgCuxH8uuKYjIr/GOKkoWkfPA9wEXgDHmGeCHwHMiUgMI8GfGmA6r6lFKzR9Op5O/+Iu/0PkKAWBZKBhjHpnh+SbgD6x6f6XU/BYVFRXoEmxJZzQrpZSaoKGglFJqgoaCUkqpCRoKSimlJmgoKKWUmqChoJRSaoKGglJKqQlijAl0DbMiIu1AY6DrWECSAZ00qIKR/tucW0uMMTOuEzTvQkHNLREpN8aUBboOpabSf5uBoaePlFJKTdBQUEopNUFDQT0b6AKUugr9txkAek1BKaXUBO0pKKWUmqChYAMiYkTkHybd/1MR+ctJ958SkaO+28cicntAClW2IV77RGT7pMe+KCJviciYiFROuv257/l7ReSwiFSJyBEReTpwf4KFS08f2YCIDAHNwE3GmA4R+VMg2hjzlyJyL/BXwD2+59YC/w7cbIxpCWDZaoETkVXAy0Ap4MC7Z/s2oMoYEz2lrQvv/KSbjTHnRSQMyDHGNPi57AVPewr24MZ70e6/TvPcnwH/9/iud8aYQ8DzwLf8V56yI2NMLfA63n+D3wd+YYw5eZXmMXg3Bev0HTusgWANy3ZeU0HnfwHVIvK3Ux4vAiqmPFYOPOGXqpTd/RVwCBgBxieqRYhI5aQ2/8MY86KI7AQaRWQP8Abwa2OMx7/lLnwaCjZhjLkoIr8A/hgYnKG5AHpeUVnOGHNJRF4E+o0xw76HB40xa6Zp+3URWQ3cBfwpcDfwpN+KtQk9fWQv/wR8DZi8+e0RYN2Udmt9jyvlDx7fbUbGmBpjzD/iDYT/ZGlVNqWhYCPGmC7gJbzBMO5vgb8RkSQAEVmD99vXj/1eoFJXISLRInLHpIfWoAtjWkJPH9nPPwD/1/gdY8xOEckAPhIRA/QBXzHGNAeqQGV7U68pvAX8NfDfROQneE9/XkJPHVlCh6QqpZSaoKePlFJKTdBQUEopNUFDQSml1AQNBaWUUhM0FJRSSk3QIalKXYVv7sYe3900YAxo990vAarw/g7VA08YYwZE5AzeYb1jgHt8j2ER+SFwP95JWm3Ak8aYJj/9UZS6bjokVanr4FtqvN8Y8/e++/3jK3mKyAtAhTHmR75QKBtfYHDS8bHGmIu+n/8YKDTGfNOffwalroeePlLqs/sQWHqtBuOB4BOFri2lgpSGglKfgYg4ge1Aje8hA7wtIhUi8tSUtn8tIueAR4Hv+bdSpa6PhoJSN2Z8KYZy4CzwM9/jtxlj1uINim+JyKbxA4wx/90YkwW8wKSlRpQKJhoKSt2YQWPMGt/tvxhjRgDGLx4bY9qAfwNunubYX6ErfKogpaGg1BwRkSgRiRn/GfgDoNZ3f9mkpvcBR/1foVIz0yGpSs2dRcC/iQh4f7d+ZYx5y/fc/xSRFXiHpDYCOvJIBSUdkqqUUmqCnj5SSik1QUNBKaXUBA0FpZRSEzQUlFJKTdBQUEopNUFDQSml1AQNBaWUUhM0FJRSSk34/wE5ekSiDmN2swAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "df_data.rename(columns={ \"data1\": gene_expre, \"data2\": gene_mutation }, inplace=True)\n", "sns.violinplot( x=df_data[gene_mutation], y=df_data[gene_expre], palette=\"Pastel1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## BigQuery to Compute statistical association \n", "\n", "The T-score (T), assuming the distributions of gene expression with and without mutation have unequal variances, is computed by using the following equation: \n", "\n", "$$T = \\frac{ \\bar{g}_y - \\bar{g}_n }{\\sqrt{ \\frac{s_y^2}{N_y} + \\frac{s_n^2}{N_n} } }$$\n", "where\n", "\n", "- $\\bar{g}_y$ and $\\bar{g}_n$ are mean gene expression of participants with and without mutation.\n", "- $N_y$ and $N_n$ are the number of participants in the group with and without mutation.\n", "- $s_y^2$ and $s_n^2$ are the variance of gene expression for the participants with and without mutation, respectively.\n", "\n", "Since the Somatic mutation table contains information of positive somatic mutation only, the averages and variances needed to compute the $T$ score are computed as a function $S_y=\\sum_{i=1}^{N_y}{g_i}$ and $Q_y=\\sum_{i=1}^{N_y}{g_i^2}$, the summs over the gene expression and squared gene expression for articipants with somatic mutation. The following query string computes $S_y$ and $Q_y$: " ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "summ_table = \"\"\"\n", "summ_table AS (\n", "SELECT \n", " Symbol,\n", " COUNT( n1.ParticipantBarcode) as Ny,\n", " SUM( n1.data ) as Sy,\n", " SUM( n1.data * n1.data ) as Qy\n", " \n", "FROM\n", " table1 AS n1\n", "INNER JOIN\n", " table2 AS n2\n", "ON\n", " n1.ParticipantBarcode = n2.ParticipantBarcode\n", "GROUP BY Symbol\n", ")\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After computing $S_y$ and $Q_y$ we can compute the mean and the variance as :\n", "$$\\bar{g}_y =\\frac{S_y}{N_y}$$ \n", "$$s_y^2 = (N_y-1)^{-1} \\left[ Q_y - \\frac{S_y^2}{ N_y} \\right]$$ \n", "\n", "To compute $S_n$ and $Q_n$, we first compute the sums of the gene expression and squeared gene expression using all the samples and then substract $S_y$ and $Q_y$. The following query uses this approach to compute the necessary variances and means, and then computes $T$ score. The $T$ score is only computed if the variances are greater than zero and if the number of participants in each group is greater than a user defined threshold." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " in runQuery ... \n", " the results for this query were previously cached \n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
NyNnavg_yavg_ntscore
02482612.4951342.89468624.426213
\n", "
" ], "text/plain": [ " Ny Nn avg_y avg_n tscore\n", "0 248 261 2.495134 2.894686 24.426213" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "query_tscore = \"\"\"\n", "SELECT \n", " Ny, Nn,\n", " avg_y, avg_n,\n", " ABS(avg_y - avg_n)/ SQRT( var_y /Ny + var_n/Nn ) as tscore\n", "FROM (\n", "SELECT Ny, \n", " Sy / Ny as avg_y,\n", " ( Qy - Sy*Sy/Ny )/(Ny - 1) as var_y, \n", " Nt - Ny as Nn,\n", " (St - Sy)/(Nt - Ny) as avg_n,\n", " (Qt - Qy - (St-Sy)*(St-Sy)/(Nt - Ny) )/(Nt - Ny -1 ) as var_n\n", "FROM summ_table as n1\n", "LEFT JOIN ( SELECT Symbol, COUNT( ParticipantBarcode ) as Nt, SUM( data ) as St, SUM( data*data ) as Qt\n", " FROM table1 GROUP BY Symbol\n", " ) as n2\n", "ON n1.Symbol = n2.Symbol \n", ")\n", "WHERE\n", " Ny > {0} AND Nn > {0} AND var_y > 0 and var_n > 0\n", "\"\"\".format(str(MinSampleSize))\n", "\n", "sql = ( sql_data + ',\\n' + summ_table + query_tscore )\n", "df_tscore = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n", "df_tscore" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To test our implementation we can use the 'ttest_ins' function available in python:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " DRG2 \n", " mean count\n", "TP53 \n", "NO 2.894686 261\n", "YES 2.495134 248\n", "Ttest_indResult(statistic=24.426213134587776, pvalue=4.38554327085175e-84)\n" ] } ], "source": [ "print( df_data.groupby(gene_mutation).agg(['mean', 'count']) )\n", " \n", "Set1 = df_data[df_data[gene_mutation]=='NO']\n", "Set2 = df_data[df_data[gene_mutation]=='YES']\n", " \n", "print( stats.ttest_ind(Set1[gene_expre], Set2[gene_expre], equal_var=False ) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.5.4" } }, "nbformat": 4, "nbformat_minor": 2 }