{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Regulome Explorer ChiSquare test for categorical features\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 how Regulome Explorer uses the Chi-square test to compute statistical associations between two categorical features. This test is used when one of the categorical features has more than two categories. Fisher's exact test is used for the special case in which both features has only two categories. The Fisher's exact test is described in another notebook.\n", "\n", "We will use clinical data and Somatic mutations for this test, both of these features are available in BigQuery Tables. Details of the Chi-sqaure ttest can be found in the following link: https://en.wikipedia.org/wiki/Chi-squared_test\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Authenticate with Google (IMPORTANT)\n", "Our first step is to authenticate with Google -- you will need to be a member of a Google Cloud Platform (GCP) project, with authorization to run BigQuery jobs in order to run this notebook. If you don't have access to a GCP project, please contact the ISB-CGC team for help (www.isb-cgc.org)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Import Python libraries" ] }, { "cell_type": "code", "execution_count": 2, "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", "import seaborn as sns\n", "import re_module.bq_functions as regulome" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Userdefined Parameters\n", "The parameters for this experiment are the cancer type, the name of the clinical feature, the name of the gene for which mutation can be extracted, and the minimun number of participant for categories to be considered. The clinical feature must be categorical. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "cancer_type = 'BRCA'\n", "clinical_name = 'histological_type'\n", "mutation_name = 'CDH1'\n", "MinSampleSize = 10\n", "\n", "bqclient = bigquery.Client()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data from BigQeury tables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clinical data from the BigQuery. The following string query will retrieve clinical data fromthe 'pancancer-atlas.Filtered.clinical_PANCAN_patient_with_followup_filtered' table available in pancancer-atlas dataset. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "query_table1 = \"\"\"table1 AS (\n", "SELECT\n", " symbol,\n", " avgdata AS data,\n", " ParticipantBarcode\n", "FROM (\n", " SELECT\n", " '{0}' AS symbol, \n", " {0} AS avgdata,\n", " bcr_patient_barcode AS ParticipantBarcode\n", " FROM `pancancer-atlas.Filtered.clinical_PANCAN_patient_with_followup_filtered`\n", " WHERE acronym = '{1}' AND {0} IS NOT NULL \n", " )\n", ")\n", "\"\"\".format(clinical_name, cancer_type)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Somatic mutation data from Bigquery table. 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": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "query_table2 = \"\"\"table2 AS (\n", "SELECT\n", " symbol,\n", " ParticipantBarcode\n", "FROM (\n", " SELECT\n", " Hugo_Symbol AS symbol, \n", " ParticipantBarcode AS ParticipantBarcode\n", " FROM `pancancer-atlas.Filtered.MC3_MAF_V5_one_per_tumor_sample`\n", " WHERE Study = '{1}' AND Hugo_Symbol = '{0}'\n", " AND FILTER = 'PASS' \n", " GROUP BY\n", " ParticipantBarcode, symbol\n", " )\n", ")\n", "\"\"\".format( mutation_name , cancer_type )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following query combines the two tables based on Participant barcodes. Data of participants for which one feature is missing are not being used. Nij is the number of participants for each pair of categories. data1 is the categorical data fo the clinical feature specified by the user, and data is binary data which is 'YES' for pariticpants with mutation in the gene especified by the user. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "query_summarize = \"\"\"summ_table AS (\n", "SELECT \n", " n1.data as data1,\n", " IF( n2.ParticipantBarcode is null, 'NO', 'YES') as data2,\n", " COUNT(*) as Nij\n", "FROM\n", " table1 AS n1\n", "LEFT JOIN\n", " table2 AS n2\n", "ON\n", " n1.ParticipantBarcode = n2.ParticipantBarcode\n", "GROUP BY\n", " data1, data2\n", ") \n", "\"\"\".format(str(MinSampleSize) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point we can take a look at output table, where the column **Nij** is the number of participants for each pair of categorical values." ] }, { "cell_type": "code", "execution_count": 7, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
data1data2Nij
0Infiltrating Carcinoma NOSNO1
1Infiltrating Ductal CarcinomaNO768
2Infiltrating Ductal CarcinomaYES9
3Infiltrating Lobular CarcinomaYES83
4Infiltrating Lobular CarcinomaNO118
5Medullary CarcinomaNO5
6Medullary CarcinomaYES1
7Metaplastic CarcinomaNO8
8Mixed Histology (please specify)NO26
9Mixed Histology (please specify)YES4
10Mucinous CarcinomaNO17
11Other specifyYES3
12Other specifyNO43
\n", "
" ], "text/plain": [ " data1 data2 Nij\n", "0 Infiltrating Carcinoma NOS NO 1\n", "1 Infiltrating Ductal Carcinoma NO 768\n", "2 Infiltrating Ductal Carcinoma YES 9\n", "3 Infiltrating Lobular Carcinoma YES 83\n", "4 Infiltrating Lobular Carcinoma NO 118\n", "5 Medullary Carcinoma NO 5\n", "6 Medullary Carcinoma YES 1\n", "7 Metaplastic Carcinoma NO 8\n", "8 Mixed Histology (please specify) NO 26\n", "9 Mixed Histology (please specify) YES 4\n", "10 Mucinous Carcinoma NO 17\n", "11 Other specify YES 3\n", "12 Other specify NO 43" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sql = ( 'WITH\\n' + query_table1 + ',' + query_table2 + ',' + query_summarize +\n", "\"\"\"SELECT * FROM summ_table \n", " ORDER BY data1\n", "\"\"\")\n", "\n", "df_results = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n", "df_results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The table shows that **data2** (Gene mutations) has two categories and **data1** (Clinical feature ) in this case has 8 categories. We can use python to visualize the populations in each category. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/borisaguilar/anaconda3/lib/python3.7/site-packages/matplotlib/tight_layout.py:176: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations. \n", " warnings.warn('Tight layout not applied. The left and right margins '\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEmCAYAAABBMrbjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xm8XdP9//HXW2KIhKSGkqY0qKFEhIRSUyj9mkpVWjRFdEC/bQ0/w5f22zRoKdFq1VR8CaqomVQNjYg5JGQk+BL9NuZ5piSf3x9rnWQ7zr33nGTfKff9fDzu4+6z9lprr7PvfZzPWWvvvZYiAjMzs7Is0d4NMDOzxYsDi5mZlcqBxczMSuXAYmZmpXJgMTOzUjmwmJlZqRxYzMysVA4sZmZWKgcWMzMrlQOLdVk77bRTAP7xT1v9dBkOLNZlvfLKK+3dBLPFkgOLmZmVyoHFzMxK5cBiZmalcmAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzEq1WAUWSe/UkWdrSTMlTZHUT9LVOX2opLGF7a8sxPEHSdql8Hp3Scc2Wk8z9R8laZakGZKmStq/pHrvK6OeRSXpGUnXFF4PkzSm8PobkqblczBd0jcK+zaXNDH/XR+TNKptW29mFd3buwHtYDhwWkRclF8Pq5FnKPAO8KkPXEndI+LjJuoeBAwBbgaIiBuBGxe1wfm4hwA7AptFxFuSegPfaKFYsXyT7Y6IhoNoKxoiaYOImFlMlLQRcBqwY0TMlrQGcLukpyNiGnAx8O2ImCqpG7Bu2zfdzGAx67FU5B7HnZKuzt9uL1PyA+DbwMic1l/SjKqy/YFDgCPyt9+tJY2R9DtJ44FTJG0m6T5Jj+Tf60paCjgB2DuX21vSCEln5nrHSDoj539a0rCcvoSks3Mvaqykmyv7qvwM+M+IeAsgIt6MiItzHSMlPZR7MudJUk6/U9JJkiYAh0laRdJ1ubcztdIrq/T0mjpved9X8/udLulCSUvn9GfyMe6XNEnSJpJulfRUDoZI6iVpnKSHc/k9mvnznZbfa7WjgJMiYnZ+/7OBk4Gj8/7PAs/nfXMj4tFmjmFmrWixDCzZxsDhwPrAmsCWEXEBqQdxdEQMr1UoIp4BzgVOj4hBEXF33rUOsENEHAnMAraJiI2BkaQPvH/n7StzuStrVN8X2ArYDfhNTvsm0B/YEPgBsEV1IUnLActFxFNNvNczI2LTiBgA9Mj1V/SJiG0j4rfAGcCEiNgI2ASYWaOuT503ScsAY4C9I2JDUk/3R4Uy/4qILYC7c75hwOakQAvwAbBnRGwCbAf8thKwavgrsImkL1albwBMrkqblNMBTgcez4Hz4NxmM2sHi3NgeTAi5kTEPGAK6cN7UVwVEXPzdm/gqtzbOZ0FH24tuT4i5uVv06vktK1y3fMi4gVgfI1yovn1HLbL1xemA9tXtacY4LYHzoH53+rfrFFXrfO2LjA7Ip7IeS4GtimUqQz3TQcmRsTbEfEy8IGkPrn9J0maBvwD6Fd4/9XmAqOB46rSa52D+WkRcQJpGPI24DvALbUql3RQ7llNevnll5togpktisU5sHxY2J7Lol9PerewfSIwPvcQvg7U++242CZV/W5SHv56V9Ka1fvyN/OzgWG5N3F+VXverS7TQBsr562lNlbKzKsqPy+XHw6sDAyOiEHAizR/zi4lBa7VC2kzSYGjaBNg/pBXRDwVEecAXwU2krRidcURcV5EDImIISuvvHILb8vMFsbiHFgWxdvAcs3s7w08m7dHNFCulnuAvfK1llVINw7UcjJwlqTlASQtL+kgFnxAvyKpF7VvRqgYRx7CktStUlcdZgH9C8NT+wET6iwL6Xy9FBEfSdoO+EJzmSPiI1JP8PBC8mnAcfkaWOVa2M+A3+bXuxaG19YmBcU3GmijmZXEgaW2m4A9Kxfva+w/FThZ0r1At0L6eGD9ysX7Oo91DTAHmAH8CZgI1BqiOifX/1AegpsAvBcRb5B6KdOB64GHmjnWYaRhs+mk6xV1DeFFxAfAgaThv+mknsi59ZTNLiPd7TWJ1HuZVUeZ/6HQy4yIKcB/ATdJmkX6Gx2T0yEFu8clTSH1eIYXhi7NrA0pokstxdwhSeoVEe/koZsHSTcavNDe7VrcDRkyJCZNmtTezbCuo8Vh78VFV3yOpSMamy9yLwWc6KBiZp2ZA0sHEBFD27sNZmZl8TUWMzMrlQOLmZmVyoHFzMxK5cBiZmalcmAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErl2Y2ty+q56hqx3n7Ht3czbDEyefT+ze3uMrMbu8diZmalcmAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErV4QOLpHfqyLO1pJmSpkjqJ+nqnD5U0tjC9lcW4viDJO1SeL27pGMbraeJusdImi1pqqQnJF0iqd9C1tVH0n/WmbfmOZW0qqQrJD0l6VFJN0taZ2HaU1VvaefMzDq+Dh9Y6jQcOC0iBkXEsxExrEaeoUDNwCKpezN1DwLmB5aIuDEifrMoja1ydERsBKwLPAKMl7TUQtTTB6grsNQiScB1wJ0RsVZErA/8DFil3vKSav4/tcI5M7MOrNMEltzjuFPS1ZJmSbosf5j9APg2MDKn9Zc0o6psf+AQ4Ijcq9k69xZ+J2k8cIqkzSTdJ+mR/Hvd/AF/ArB3Lre3pBGSzsz1jpF0Rs7/tKRhOX0JSWfnXtTY/M2/VrCbL5LTgReAnXM983sWkoZJGpO3V5F0Xe7pTM09sd8Aa+V2jpbUS9I4SQ9Lmi5pjxZO8XbARxFxbqFNUyLi7qbqyuf6MUlnAw8Dq0naKeebKmlczlfPOVNu94x8jL0Lf/cJkv6ae3W/kTRc0oM531o539clTcx/v39Iqisgmln5mvum3hFtDGwAPAfcC2wZERdI2goYGxFX5yDyCRHxjKRzgXci4jQASd8H1gF2iIi5kpYHtomIjyXtAJwUEXtJGgkMiYif5HIjqqrvC2wFrAfcCFwNfBPoD2wIfBZ4DLiwzvf4cK7rhmbynAFMiIg9JXUDegHHAgMiYlBuZ3dgz4h4S9JKwAOSboymJ4cbAExuYt8HterK+9YFDoyI/5S0MnA+6TzOlrRCE/U1dc4GARsBKwEPSbor598I+BLwGvA0cEFEbCbpMOCnwOHAPcDmERH5y8YxwJHVB5Z0EHAQQL/eS3LdcqNZfeT0JpppZgujswWWByNiDoCkKaQP73sWob6rImJu3u4NXCxpbSCAJeus4/qImAc8WviWvFWuex7wQu4V1aueieq2B/YHyO1/U9JnatRzkqRtgHlAP9Kw1gsNtKWlugD+GREP5O3NgbsiYnZu22tN1NfUObs8v58XJU0ANgXeAh6KiOcBJD0F3JbLTCf1tAA+D1wpqS+wFDC71oEj4jzgPICB/Xp4BlazVtBphsKyDwvbc1n0wPhuYftEYHxEDAC+DiyzEG1S1e+FsTGphwMpwFXU256K4cDKwODci3mxhTpmAoMXoq7iOVRVm5vS6Dkr5p9XeD2PBf8DfwTOjIgNgYNp/HyZWUk6W2BZFG8DyzWzvzfwbN4e0UC5Wu4B9srXWlYh3TjQrHyN4VDSMNEtOflFSV9Suii+ZyH7OOBHuVy3PIxX3c7ewEsR8ZGk7YAvtNCEO4ClJf2w0KZNJW3bQF33A9tKWiOXb2oorJa7SNeyuuUhtW2ABxsoX/z7HdBAOTMrWVcKLDcBe+aL21vX2H8qcLKke4FuhfTxwPqVi/d1HusaYA4wA/gTMBF4s4m8oyVNBZ4gDf1sFxH/zvuOBcaSPvSfL5Q5DNhO0nTSdZENIuJV4N588Xs0cBkwRNIkUo9jVnMNztde9gR2VLrdeCYwinQ9q666IuJl0vWLa/N7urK5Y1a5DpgGTM3v95iIaGTYbhRwlaS7gVcaKGdmJfNCX61EUq+IeEfSiqRv3ls2+EFprWxgvx4x9uAv+uK9tZUus9BXZ7t435mMldSHdCH5RAcVM+sqHFhaSUQMbe82mJm1h650jcXMzNqAA4uZmZXKgcXMzErlwGJmZqVyYDEzs1I5sFiXtVTfDfwMi1krcGAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVZze2LqvnqmvEevsd3yp1Tx69f6vUa51al5nd2D0WMzMrlQOLmZmVyoHFzMxK5cBiZmalcmAxM7NSObCYmVmpHFjMzKxUHS6wSHqnjjxbS5opaYqkfpKuzulDJY0tbH9lIY4/SNIuhde7Szq20XqaqHuMpGF15p3/Xhqof5SkoxaudfPr2EzSXZIelzRL0gWSll2UOnO9J0jaYVHrMbOOr3t7N2AhDQdOi4iL8utaH9ZDgXeA+6p3SOoeER83UfcgYAhwM0BE3AjcuKgN7oiqz4OkVYCrgH0i4n5JAvYClgPea7S+oogYWVKzzayD63A9lor8jf1OSVfnb86XKfkB8G1gZE7rL2lGVdn+wCHAEblXs3XuLfxO0njglPzN/D5Jj+Tf60paCjgB2DuX21vSCEln5nrHSDoj53+60vuQtISks3MvaqykmxvomUjSaEkzJE2XtHdh9/KSrpP0qKRzJS2Ry7xTKD9M0pga9f5Q0kOSpkq6ptLrqD4PVcV+DFwcEfcDRHJ1RLxY63zl+kZIukrSTcBtOe2Y/F6mSvpN4biV8/WMpOMlPZzzrZfTV5B0vaRpkh6QNDCnj5J0saTbctlvSjo1l71F0pI538j8nmdIOi8HRjNrYx02sGQbA4cD6wNrAltGxAWkHsTRETG8VqGIeAY4Fzg9IgZFxN151zrADhFxJDAL2CYiNgZGAidFxL/z9pW53JU1qu8LbAXsBvwmp30T6A9sCPwA2KKB9/hNUi9pI2AHYLSkvnnfZsCRud61ct56XRsRm0bERsBjwPcL+4rnoWgAMLmJ+j51vgr7tgAOiIjtJe0MfAP4cj72qU3U90pEbAKcA1SG744HHomIgcDPgEsK+dcCdgX2AP4MjI+IDYH3czrAmfk9DwB6kP5GZtbGOvpQ2IMRMQdA0hTSh/c9i1DfVRExN2/3Bi6WtDYQwJJ11nF9RMwDHs1DR5ACzVU5/YXcG6jXVsDluV0vSpoAbAq8RXr/TwNIujznvbrOegdI+hXQB+gF3FrYVzwP9WrufN0eEa/l7R2AiyLiPYBCerVr8+/JLAiYW5GG3oiIOyStKKl33vf3iPhI0nSgG3BLTp9O+r8A2E7SMcCywArATOCm4kElHQQcBNCv95Jct9zoOt/+p3n1SbPaOnqP5cPC9lwWPRC+W9g+kfStdwDwdWCZhWiTqn4vjObKVs8QGjXSm2r3GOAn+Vv98VX53q1ZIn0QD25iX3Pnq1ifarS7lsp5LP5da52LSl0fAuTg/VEsmD11HtBd0jLA2cCw/J7Pp8a5iYjzImJIRAxZoWe3OpppZo3q6IFlUbxNuujclN7As3l7RAPlarkH2Ctfa1mFdONAve4iXdPpJmllYBvgwbxvM0lr5Gsre7Ogt/aipC/l9D2bqHc54Pl8/aHmkGENZwIHSPpyJUHSdyWtStPnq9ptwPcK13RWqPPYkM7F8FxuKGm47K06y1aCyCuSelH7hg4zawOLc2C5CdizcvG+xv5TgZMl3UsaWqkYD6xfuXhf57GuAeYAM4A/AROBN5vI+ydJc/LP/cB1wDRgKnAHcExEvJDz3k+6jjMDmJ3zAhwLjM35n2/iOL/I7biddH2kRRHxIrAPcJrS7caPAVuThuWaOl/VddxCugY2KQ9fNnL78yhgiKRppPd9QL0FI+INUi9lOnA98FADxzWzEnk9lpJI6hUR70hakdTj2LIQIKwDGtivR4w9+IsLXd7XWKxBXeYuxY5+8b4zGSupD7AUcKKDipl1VXUHlvxMwHBgzYg4QdLqwKoR8WALRbuEiBja3m0wM+sIGrnGcjbpeYV98+u3gbNKb5GZmXVqjQyFfTkiNpH0CEBEvK70pLqZmdl8jfRYPpLUjfxcQb41dl6rtMrMzDqtRgLLGaTbXVeR9GvSMxUnNV/EzMy6mrqHwiLiMkmTga/mpG9ExGOt0ywzM+usGr3deFnSw3FBmuTPrNNaqu8GrD5yUns3w2yxU/dQmKSRwMWkyf1WAi6S9N+t1TAzM+ucGumx7AtsHBEfAOR1Nh4GftUaDTMzs86pkYv3z/DJ2WKXBp4qtTVmZtbpNdJj+RCYKel20jWWHYF7JJ0BEBGHtkL7zMysk2kksFzHgtl1Ae4stylmZrY4aCSwvArcnBdaMuv0Hpvzans3wWyx1Mg1ln2AJyWdKulLrdUgMzPr3OoOLBHxXWBj0gX7iyTdL+kgSY2utmhmZouxhlaQzMvEXgNcAfQlLYv7sKSftkLbzMysE2rkAcndJV1HWg53SWCziNgZ2IjGlp81M7PFWCMX74cBp0fEXcXEiHhP0vfKbZaZmXVWjQyFPV8dVCSdAhAR40ptlZmZdVqNBJYda6TtXFZDzMxs8dDiUJikHwH/CawlaVph13LAva3VMDMz65zq6bH8Bfg6cEP+XfkZnG9BBkDSZ1qlhR2YpJB0aeF1d0kvSxrbYD13ShrSQp5Rko7K22MkDVu4VrfYlv0lzZA0U9KjlWOWUO/NkvqUUZeZdWwt9lgi4k3gTdLsxs0ZB2xSRqM6kXeBAZJ6RMT7pOHCZ9u5TZ8iqXtEfFxHvp2Bw4GvRcRzkpYB9ivjOBGxS90NNrNOraHnWFqgEuvqTP4O7Jq39wUur+yQ1FPShZIekvSIpD1yeg9JV0iaJulKCoumSXqnsD1M0pjmDi5pZK5/hqTzJCmn3ynpJEkTgJ9Lmi1pybxveUnPVF4XHAccFRHPAUTEBxFxfi7zw3ycqZKukbRsTh8j6XeSxgOnSOol6SJJ0/P72yvne0bSSpL6S3pM0vm5V3SbpB45zyBJD+Ry11V6wfm9nC7prlx2U0nXSnpS0vxlGyRdL2lyrveguv56Zla6MgNLlFhXZ3IFsE/+dj8QmFjY93PgjojYFNgOGC2pJ/Aj4L2IGAj8Ghi8CMc/MyI2jYgBpAC1W2Ffn4jYNiKOJ00aWgmA+wDXRMRHVXUNACY3cZxr83E2Ah4Dvl/Ytw6wQ0QcCfwCeDMiNszv744ada0NnBURGwBvAHvl9EuA/8rlpgO/LJT5d0RsA5xLGpb9cW7vCEkr5jzfi4jBwBDg0EK6mbWhMgNLlxQR04D+pN7KzVW7vwYcK2kK6YN9GWB1YBvgz4Xy01h420maKGk6sD2wQWHflYXtC4AD8/aBwEUNHmeApLvzcYZXHeeqiJibt3cAzqrsiIjXa9Q1OyKm5O3JQH9JvUmBcEJOv5h0nipuzL+nAzMj4vmI+BB4Glgt7ztU0lTggZy2dvWB8zREkyRNWmlpz6dq1hoaXfO+OV11KAzSh95pwFCg+C1ZwF4R8Xgxcx6taqqHV0xfpok8lXqWAc4GhkTEvySNqirz7vxKI+7Nw1DbAt0iYkaNKmeSek+1ehljgG9ExFRJI0jv9VPHIb3nlnqvHxa251IYCqyjzLyq8vOA7pKGkoLaFvmh3Tupcf4i4jzgPIAhQ4Z01V62WatqscciaYXmfgpZv9qK7ezoLgROiIjpVem3Aj8tXPfYOKffRfrWj6QBpCG0ihclfUnSEqS52JpT+eB8RVIv0uwIzbmEdA2oqd7KycCpklbNbVtaUmUBt+WA5/N1meHNHOM24CeVF6rzbsF8k8jrkrbOSfsBE5opUq038HoOKusBmzdQ1sxKVE+PZTLpG2itHkkAawJExGsltqtTiYg5wB9q7DoR+D0wLQeXZ0jXQM4hzRA9DZgCPFgocywwFvgXMAPo1cxx35B0Pml46BngoRaaehnwKwo3GFTVd7OkVYB/5PYGKWhCunYyEfhnPl5Ts1r/CjhL0gxSb+R44NoW2lVxAHBuvjHgaRYM3dXjFuCQfE4fJw2HmVk7UIRHA7oKpWdf9oiIum8hXpwNGTIkJk2a1N7NsK6jy1wuaOgaSx7WWJvC2HX1/GHWMUn6I2kKHj9PYmatqu7AIukHwGHA50nDN5sD95PuRLIOLiK8Zo6ZtYlGbjc+DNgU+GdEbEdaTfLlVmmVmZl1Wo0Elg8i4gNIdwtFxCxg3dZplpmZdVaNBJY5SpMIXg/cLukG4LnWaZaZmVVIWjVPA/WU0uSwN0taR9L7StNFPSbpQUkHFMqMkHRmVT3zJ7yV9GtJ/1JhGqmy1H2NJSIqz1SMyvNC9Sbd4mlmZq0k3/p/HXBxROyT0wYBqwBPRcTGOW1N4FpJS0REPTNr3AScCTxZdpsbWfN+c0nLAeRpN8aTrrOYmVnr2Q74KCLOrSTkKZH+VcwUEU8D/w84lDpExAMR8XyZDa1o5Hbjc/jktPjv1kgzM7NyNTc5bLWHgfUKr/eWtFXh9RdLa1UzGgksisLTlBExT1KZc42ZmdmiqX4I88qIKE6xdGdbNKKRi/dPSzpU0pL55zDStBtmZtZ6KpPD1mNj0rIW7aqRwHII8BXSColzgC8DXkzJzKx13QEsLemHlQRJmwJfKGaS1J80y/of27JxtTRyV9hLpAWizMysjURESNoT+L2kY4EPSJPOHg6sJekR0jRbbwN/rPOOMCSdCnwHWFbSHOCCiBhVRptbnIRS0jERcWqea+pTmSOirjsQzDoaT0JpbazLTEJZz1BYZbxuEunOhOofs07psTmvMvjoS9q7GWaLnRaHwiLipvz74tZvjpmZdXaNzG58E58eCnuT1JP5U2UeMTMz69oaut0YeAc4P/+8BbwIrJNfm5mZNfSA5MYRsU3h9U2S7oqIbSTNLLthZmbWOTXSY1lZ0uqVF3l7pfzy36W2yszMOq1GAsuRwD2SxudpAe4GjpbUE/CFfTOzDk5SSPpt4fVRkkYVXh8kaVb+ebBqnrG6NfKA5M2S1iZNcCZgVuGC/e8X5uBmZl3V4KMvaf4hwgZNHr1/Pc/JfAh8U9LJEfFKcYek3YCDga0i4hVJmwDXS9osIl5opC2NTJu/ZD7oL4D/Bn6Q08zMrHP4GDgPOKLGvv8Cjq4EnIh4mDQa9eNGD9LIUNg5pInQzs4/g3PaYit3Gy8tvO4u6WVJY1soN0jSLot47IVa1U3S4ZKWLby+Oa/8WW/5nSVNyivSzZJ02sK0o0a9F0hav4y6zGyRnAUMl9S7Kn0DPv3Q+6Sc3pBG7grbNCI2Kry+Q9LURg/YybwLDJDUIyLeB3YkTcLZkkHAEODm1mxcEw4H/gy8BxARdQc4SQNIK8rtGhGz8rIIdU80Kql7RHxca19E/KDeesys9UTEW5IuIS0I9n4L2UWNqbxa0kiPZa6kteYfLS2DObfRA3ZCfwd2zdv7ApdXdkjqKelCSQ/ldaf3kLQUcAJpgZ0pkvaWtJmk+3Ke+yStm8uPkHSDpFskPS7pl9UHl9RL0jhJD0uaLmmPwrH/JmmqpBn5OIcCnwPG5+WjkfSMpJXy9v6SpuUyl1YfCzgG+HVEzAKIiI8j4uxc9uuSJub38A9Jq+T0UZLOk3QbcImkbpJOy22dJumnOV9xre13lNbbnirpgUJdX8jvdVr+vXpOHyPpnHzjyNOSts3n/TFJYwrn6pzc25op6fiF+FubdRW/B74P9CykPcqnp+ffJKc3pJHAcjTpA+tOSRNIUzkf2egBO6ErgH0kLQMMBCYW9v0cuCMiNiUtHzoaWBIYSVpgZ1BEXAnMArbJa1OPBE4q1LEZMJzUy/lW5cO34ANgz4jYJB/jt5IE7AQ8FxEbRcQA4JaIOAN4DtguIrYrViJpg9ze7XPP87Aa77W5leruATbP7+EKUhCqGAzsERHfIfVw1iA99zQQuKxGXT2BB3I77gIq04GfCVxSKHdGocxngO1JY8M3AaeTuugbKq3/DfDziBhC+jttK2lg9YHzXS+TJE1aael5TB69fxNv12zxFRGvAX8lBZeKU4FTJK0IaUgfGEG69NGQRu4KG5fvCluXBXeFfdjoATubiJimtM7Bvnx6aOtrwO6SjsqvlwFW59N6Axfn8xek4FNxe0S8CiDpWmAr0rhmhYCTJG0DzAP6AasA04HTJJ0CjI2Iu1t4K9sDVxcuzL3WQv5qnweulNQXWAqYXdh3Yx4qBNgBOLcyJNbEcf4NVK5TTSYNMQJsAXwzb19K+kevuClPHz4deDEipgMoPZzbH5gCfFvSQaT/677A+sC04oEj4jzSxUuGDBlS6l05Zp3Mb4H5q0tGxI2S+gH3SQrSNPzfjYjnG624xcAi6ZtN7FpLEhFxbaMH7YRuJC2gMxRYsZAuYK+IeLyYWdKXq8qfCIyPiD1zkLqzsK/6w6369XBgZWBwRHwk6RlgmYh4QtJgYBfgZEm3RcQJzbyHesZKKyvV1bp29kfgd/mfbygwqrDv3QaP81Fhmeu5NP1/WKyn8iVmXmG78rq7pDWAo0jXAl/PQ2TLtNAOs3ZT5+3BpYqIXoXtF4Flq/afQwk3ZdUzFPb1Zn52W9QGdBIXAidUviUX3Ar8NA9NIWnjnP42sFwhX28WXPQfUVXHjpJWkNQD+AZwb9X+3sBLOahsR141TtLngPci4s+koLdJE8euGEf6Rl/p5q5QI89o4GeS1sl5lpD0/2q8hwNqlK24DTgkX/hv6jhNuY8Fi8kNJw2/1Wt5UoB7M1+z2bmBsmZWonqmzT+wLRrSkUXEHOAPNXadSLoINi0Hl2dIwXY8cKykKcDJpCGdi/OH9B1VddxDGvb5IvCXiKheeeoy0rxsk0jDPbNy+obAaEnzgI+AH+X084C/S3q+eJ0lImZK+jUwQdJc4BGqglwe9jscuFzpluUA/pZ3jwKukvQs8ADpOkotF5AmJp0m6SPSBKVnNpG32qHAhZKOBl4G6v7fi4ipSivpzSRNmFodoM2sjbS4guT8jOme518ClYkoJ5C+xb/ZSm1b7EkaAQyJiJ+0lNfK5xUkrY15BckaLiQNs3w7/7wF1LW2spmZdR2NPCC5VkTsVXh9fB7qsYUUEWOAMe3cDDOzUjXSY3lfhZkuJW1Jy09tmplZF9NIYDlQVJZnAAAWRklEQVQEOCs/yf0M6YLswa3SKjMzK5WSeyTtXEj7ttLMH3OVZgqp/Byb9++WZ9uYKulRSXV95jcyFPZWRGwkaXmYP99MU3cGmZlZM/7vhA1LfUB39ZHTm705ID9gfAjp7s7xQDfg16RZPKZGxKBifqXZ688DNouIOZKWJj2M3KJGAss1wCYR8VYh7Wo+PbeMmZl1QBExQ9JNpCnye5KmUHoqP4pXbTlSjHg1l/0QeLxWxmr1PHm/HmlOpt5VT+Evj59sNjPrbI4HHiZNrVSZm7BH1c1YJ0fElZJuBP4paRxpGqbLI2JeSweop8eyLumhvz6kp+0r3mbB5IFmZtYJRMS7kq4E3inM9/h+9VBYzvsDSRuS5gA8ijSv34iWjlHPk/c3ADdI2iIi7m/kDZiZWYc0L/+0KE9lNV1pqY3Z1BFYGrkrbE9Jy0taUmmtjFckfbeB8mZm1kkorQU1tJA0CPhnPWUbCSxfyxfudwPmkOaDOrqB8mZm1jH1qLrd+DekKWiOUVqEcArp2syIeipr5K6wyhoiu5Au4LzWxJ0EZmbWgpZuD25NETGq6nW3JrLWvbR5USOB5SZJs0hP2/+npJVJqxuamZnNV/fsxgCSPkN6UHJunlZ9+Yh4odVaZ9aKeq66Rqy33/GfSPNSxdaKuswQTz3PsWwfEXcUn2GpGgLrCitImplZneoZCtuWtDhV5RmWShensgStA4uZmc1Xz3Msv8ybPwL2Is0VUylX6lw3ZmbW+TVy8f564A3SVACVi/YOLGZm9gmNBJbPR8ROrdYSMzNbLDTygOR9ec4YMzOzJtVzV9h00pBXd+BASU8DH5Iv3kfEwNZtopmZdSb19Fh2I90RtjPwReBr+XUlvUmSIk9cVnndXdLLksbm17tXVipbVJLeqSdd0ghJZ+btQyQ1+eCCpKGSvlLHsUdJOqrRNjdC0saSLmghT39JM1qzHW2t+D8iaWVJE/OKdls3kX8pSXdJamSY18xKVM9dYXVNOtaEd4EBknpExPukKZefLdR9I3DjItS/SCLi3BayDAXeAe5r/da06GfAr9q7EW2t6n/kq8CsiDigmfz/zmtH7A1c1gZNNLMqjVxjWVh/B3bN2/sCl1d2VPUebqj0HiQdLOmyvL1WXpN5sqS788JjSFpD0v2SHpJ04sI0rNjTkHRoXtN5mqQrJPUHDgGOyJOybS3pC3lm52n59+o16hwk6YGc57o8WwGSNs1p90saXelZ5Pc0qFD+XkkDq+pcDhgYEVML7b5U0h2SnpT0qXVxJHXLx3koH/fgnN4rt/1hSdMl7ZHTe0r6m9La1jMk7Z3TB0uakM//rZL61jjWt3KZqZLuymkj8t/0ljyJ3S8L+b8r6cF8Xv8kqVtO3ym3a2oODvP/R/I5OhXYJZf7saTTC3X+UNLv8svrgeHN/OnNrBW1xXDBFcDIPPw1ELgQqDWMcRBwr6TZwJHA5jn9POCQiHhS0peBs4HtgT8A50TEJZJ+3Mzxq1dGW4HavaRjgTUi4kNJfSLiDUnnkhbDOQ1AaUnPSyLiYknfA84AvlFVzyXATyNigqQTgF8ChwMXAQdFxH1KM4dWXECaMfRwSesAS0fEtKo6hwDVQ1wD8znqCTwi6W9V+78PvBkRmyqtVX2vpNuAfwF7RsRbklYCHlBaJW4n4LmI2DW/195Ka17/EdgjIl7OwebXwPeqjjUS+I+IeFZSn0L6ZsAA4D3godzGd0m9iS0j4iNJZwPDJf0dOB/YJiJmS1qheICImCJpJDAkIn4iqScwTdIxEfERcCBwcM4+A9gUM2sXrR5YImJa/va/L3BzM/lezB8c40kffK9J6gV8BbhKC6aRWTr/3pL0wCbApcApTVT9iZXRJI1gwXKcRdOAyyRdT/rGW8sWQGVqm0tJ36Dnk9Qb6BMRE3LSxbntfYDlIqIypPYX0jUqgKuAX0g6mvSBPabGcfsCL1el3ZCHF9+XNJ70IV4MoF8DBkoall/3BtYmLXlwkqRtSAv99ANWAaYDp0k6BRgbEXdLGkAKDLfn898NeL5G++4Fxkj6K5+cieH2iHg1n5trga2Aj4HBpEAD0AN4iRQk74qI2QAR8VqN48yXV8G7A9hN0mPAknlBIvJcdv+WtFxEvF0sJ+kg0pcY+vVe0nODmbWCtrrAeSNwGumaxYrN5NsQeBX4XH69BPBGrSUzszIf0NwV2AbYnfRBv0EdZeo9fpOTz0XEe5JuB/YAvk3toPc+sEwLx65+LVLP6dZPJKbAujIwOPcYngGWiYgnJA0mTZN9cu7dXAfMjIgtmntzEXFI7k3uCkwpDO3VaqOAiyPiuKp27V4jf0suIF17mkXqERYtTY3ZtyPiPFIvmIH9evgBX7NW0BbXWCANf51Q+UZZi6TNSHeebQwcJWmNvLDYbEnfynkkaaNc5F5gn7y9SOPpkpYAVouI8cAxQB+gF/A2sFwh631Vx7ynWE9EvAm8rgV3LO0HTIiI14G3JVWG9/bhky4gDas91MQ39cdId+QV7SFpGUkrkgL2Q1X7bwV+lIezkLROHj7qDbyUg8p2wBfy/s8B70XEn0lfAjYBHgdWlrRFzrNkrYAraa2ImBgRI4FXgNXyrh0lrSCpB2nI8F5gHDBM0mdz2RUkfQG4H9hW0hqV9Brn4RMiYmI+1nf45LW7FYGX8xCZmbWxNumxRMQc0jWRmvI1gPOBAyPiOUlHAhdK2p70AX6OpP8mLTZ2BTAVOAz4i6TDgGsWsYndgD/noSwBp+drLDcBVytd4P4pcGhu19GkoakDa9R1AHCu0rICTxfyfB84X9K7wJ3Am5UCETFZ0lt8+lt3Zf+sfM2jOLTzIPA3YHXgxHze+heKXUCa1+1hpTGnl0kf7peR1taZRBo6m5XzbwiMljQP+Aj4Ub7DahhwRj433YHfAzOrmjha0tr53I0j/X0GkQLvpaSg+JeImASQ/5a35YD+EfDjiHggD1Ndm9NfIt1F2JK/AoNy8K7YjmaGXc2sdTW0HostPEm9IuKdvH0s0DciDsuvP0cKNutFxLwmyh8BvB0RF0gaReGmgo6oci0rIn7SyscZS/oiMK6Qdi1wXEQ83lzZgf16xLRn32/N5pkVdZn1WNpqKMxgV6XbZGeQ7or7FYDSLdYTgZ83FVSyc0gzHhggqY+kJ0g3ZxSDylLA9S0FFTNrPe6xWJflHou1MfdYzMzMFoYDi5mZlcqBxczMSuXAYmZmpXJgMTOzUjmwWJe1VN96Zu0xs0Y5sJiZWakcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBxczMSuVJKK3L6rnqGrHefscDeIliawuehNLMzGxhOLCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErlwFICSSHp0sLr7pJeljR2Ies7RFK7PFghaX9JMyTNlPSopKNKqvdmSX3KqMvMOrbu7d2AxcS7wABJPSLifWBH4NmFrSwizi2tZQ2QtDNwOPC1iHhO0jLAfg2U7x4RH9faFxG7lNRMM+vg3GMpz9+BXfP2vsDllR2SRhW/+eceQf+8vb+kaZKmVno9xfyS7pR0iqQHJT0haeucvoykiyRNl/SIpO1y+ghJZxaONVbSUEndJI3Jx54u6Yga7+E44KiIeA4gIj6IiPNzPT+U9FBu5zWSls3pYyT9TtJ44BRJvQrtmiZpr5zvGUkrSeov6TFJ5+de0W2SeuQ8gyQ9kMtdJ+kzhXNwuqS7ctlNJV0r6UlJvyq81+slTc71HrSwf0gzWzQOLOW5Atgnf8sfCExsqYCkDYCfA9tHxEbAYU1k7R4Rm5F6E7/MaT8GiIgNSYHs4nzspgwC+kXEgFzmohp5BgCTmyh/bURsmtv5GPD9wr51gB0i4kjgF8CbEbFhRAwE7qhR19rAWRGxAfAGsFdOvwT4r1xueuG9Avw7IrYBzgVuyO9/ADBC0oo5z/ciYjAwBDi0kG5mbchDYSWJiGm5F7IvcHOdxbYHro6IV3IdrzWR79r8ezLQP29vBfwxl5sl6Z+kD/imPA2sKemPwN+A2+psY8WA3DvoA/QCbi3suyoi5ubtHYB9Kjsi4vUadc2OiCl5ezLQX1JvoE9ETMjpFwNXFcrcmH9PB2ZGxPMAkp4GVgNeJQWTPXO+1UgB7NXigXNP5iCA1Vdf3XOEmbUC91jKdSNwGoVhsOxjPnmuKz0LAfXMAvph/j2XBV8GmprQruax8gf8RsCdpG/7F9QoOxMY3ES9Y4Cf5N7O8Sx4D5CuMVXU854+LGwX31M9ZeZVlZ8HdJc0lBTUtsi9qkeq2ghARJwXEUMiYsjKK69cx2HNrFEOLOW6EDghIqZXpT8DbAIgaRNgjZw+Dvh2ZchG0goNHOsuYHgutw6wOvB4PtYgSUtIWg3YLOdZCVgiIq4hDVdtUqPOk4FTJa2ayywt6dC8bzngeUlLVo7bhNuAn1ReVK6TtCQi3gRer1xDIt00MKGZItV6A69HxHuS1gM2b6CsmZXIQ2Eliog5wB9q7LoG2F/SFOAh4Imcf6akXwMTJM0lfcseUefhzgbOlTSd1EsZEREfSroXmE0aMpoBPJzz9wMuklT5MnFcjfbfLGkV4B+SKj2PC/PuX5CuG/0z171cE+36FXCWpBmk3sjxLBjKa8kB+T0tSxq6O7DOcgC3AIdImkYKsA80UNbMSuT1WKzLGjJkSEyaNKm9m2Fdh9djMTMzWxgOLGZmVioHFjMzK5UDi5mZlcqBxczMSuXAYmZmpXJgMTOzUjmwmJlZqRxYzMysVA4sZmZWKgcWMzMrlQOLmZmVyoHFzMxK5dmNrcvqueoasd5+xze536tLWsk8u7GZmdnCcGAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErlwLIYkfR5STdIelLSU5L+IGmpvG+QpF0KeUdJOqr9WrtwJO0u6di8vbKkiZIekbR1e7fNzBIHlsWEJAHXAtdHxNrAOkAv4Nc5yyBglyaKL8zxupVVVyMi4saI+E1++VVgVkRsHBF3t0d7zOzTHFgWH9sDH0TERQARMRc4AviepOWBE4C9JU2RtHcus76kOyU9LenQSkWSvivpwZz3T5UgIukdSSdImghs0VRDJH1L0gxJUyXdldNG5N7ULZIel/TLOo63k6SHcz3jCvWcKWkQcCqwSy73Y0mnF+r8oaTflXBezaxB3du7AVaaDYDJxYSIeEvS/wH9gZHAkIj4CaShMGA9YDtgOeBxSecAXwT2BraMiI8knQ0MBy4BegIzImJkC20ZCfxHRDwrqU8hfTNgAPAe8JCkvwHv1jqepL8D5wPbRMRsSStUvbcpkua/J0k9gWmSjomIj4ADgYOrGybpIOAggH69l/R8YGatwIFl8SGg1oyiTaUD/C0iPgQ+lPQSsAppeGkw6YMfoAfwUs4/F7imjrbcC4yR9FfS8FzF7RHxKoCka4GtgI+bON7mwF0RMRsgIl5r7oAR8a6kO4DdJD0GLBkR02vkOw84D2Bgvx6egdWsFTiwLD5mAnsVE/IQ2GrAU6QP72ofFrbnkv4fBFwcEcfVyP9BHmJrVkQcIunLwK7AlDxsBZ8OcNHU8STtXiN/Sy4AfgbMAi5qsKyZlcTXWBYf44BlJe0P8y+u/xYYExHvAW+ThrzqqWeYpM/melaQ9IVGGiJprYiYmIfMXiEFN4Adc309gG+QejZNHe9+YFtJa1TSWzpuREzMx/oOcHkjbTaz8jiwLCYiLayzJ/AtSU8CTwAfkL7BA4wnXawvXryvVc+jwH8Dt0maBtwO9G2wOaMlTZc0A7gLmJrT7wEuBaYA10TEpKaOFxEvk66FXCtpKnBlncf+K3BvRLzeYJvNrCRe6MvahKQRFG4eaMXjjAVOj4hxLeUd2K9HTHv2/dZsjlmRF/oy60wk9ZH0BPB+PUHFzFqPeyzWZbnHYm3MPRYzM7OF4cBiZmalcmAxM7NSObCYmVmpHFjMzKxUDizWZS3Vd4P2boLZYsmBxczMSuXAYmZmpfIDktZlSXobeLy929GElUgTeHY0HbVd0PHbNisidmrvhrQFBxbrsiRNiogh7d2OWjpq2zpqu8Bt60g8FGZmZqVyYDEzs1I5sFhXdl57N6AZHbVtHbVd4LZ1GL7GYmZmpXKPxczMSuXAYl2OpJ0kPS7pfyUd2w7HX03SeEmPSZop6bCcvoKk2yU9mX9/JqdL0hm5vdMkbdIGbewm6ZG8IieS1pA0MbftSklL5fSl8+v/zfv7t3K7+ki6WtKsfP626AjnTdIR+W85Q9LlkpbpKOesPTiwWJciqRtwFrAzsD6wr6T127gZHwNHRsSXgM2BH+c2HAuMi4i1gXH5Nbmta+efg4Bz2qCNhwGPFV6fQlryeW3gdeD7Of37wOsR8UXg9JyvNf0BuCUi1gM2ym1s1/MmqR9wKGnp7QFAN2AfOs45a3sR4R//dJkfYAvg1sLr44Dj2rlNNwA7kh7W7JvT+gKP5+0/AfsW8s/P10rt+TzpA3p7YCxp5cNXgO7V5xC4Fdgib3fP+dRK7VoemF1df3ufN6Af8C9ghXwOxgL/0RHOWXv9uMdiXU3lQ6BiTk5rF3kYZGNgIrBKRDwPkH9/Nmdr6zb/HjgGmJdfrwi8EREf1zj+/Lbl/W/m/K1hTeBl4KI8THeBpJ6083mLiGeB04D/A54nnYPJdIxz1i4cWKyrqbXueLvcGimpF3ANcHhEvNVc1hpprdJmSbsBL0XE5DqP35bnszuwCXBORGwMvMuCYa9a2qRt+ZrOHsAawOeAnqRhuKaO3WH+B1uLA4t1NXOA1QqvPw8819aNkLQkKahcFhHX5uQXJfXN+/sCL+X0tmzzlsDukp4BriANh/0e6COpe43jz29b3t8beK2V2jYHmBMRE/Prq0mBpr3P2w7A7Ih4OSI+Aq4FvkLHOGftwoHFupqHgLXzHTtLkS6y3tiWDZAk4H+AxyLid4VdNwIH5O0DSNdeKun757ucNgferAz9lC0ijouIz0dEf9K5uSMihgPjgWFNtK3S5mE5f6t8+46IF4B/SVo3J30VeJT2P2//B2wuadn8t620q93PWbtp74s8/vFPW/8AuwBPAE8BP2+H429FGvqYBkzJP7uQxtnHAU/m3yvk/CLdyfYUMJ1091FbtHMoMDZvrwk8CPwvcBWwdE5fJr/+37x/zVZu0yBgUj531wOf6QjnDTgemAXMAC4Flu4o56w9fvzkvZmZlcpDYWZmVioHFjMzK5UDi5mZlcqBxczMSuXAYmZmpXJgMbOaJIWk3xZeHyVpVN4+RNL+efsESTu0UzOtA/LtxmZWk6QPSHNfbRoRr0g6CugVEaPat2XW0bnHYmZN+Zi0pO4R1TskjcqBBkljJA2rzmNdlwOLmTXnLGC4pN7t3RDrPBxYzKxJkWZdvoS0kJVZXRxYzKwlvyetetizvRtinYMDi5k1KyJeA/7KgqV1zZrlwGJm9fgtsFIz+317qc3XveUsZtYVRUSvwvaLwLKF16MKWVdkMVuoyhaNeyxmttAkXUgKOPe0d1us4/ADkmZmVir3WMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBxczMSuXAYmZmpfr/zVKFFjwLfh8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "df_results.rename(columns={ \"data1\": clinical_name, \"data2\": mutation_name }, inplace=True)\n", "sns.catplot(y=clinical_name, x=\"Nij\",hue=mutation_name,data=df_results, kind=\"bar\",height=4, aspect=.7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the statistics \n", "After sumarizing the data in the table above, we are in the position to perform the Chi-Square test. Before the description of the test we need the following definitions:\n", "\n", "- **Nij** : Number of participants for each pair of categories\n", "- **N** : Total number of participants ( N = sum_ij Nij )\n", "- **Ni** : The total number of participants with category of data1(Clinical data) equal to i\n", "- **Nj** : The total number of participants with category of data2(Somatic mutation) equal to j\n", "- **E_nij** : Expected number of participants for each pair of categories under the null hypothesis. E_nij = (Ni\\*Nj)/N\n", "- **I, J** : The number of categories in data1 and data2 respectively\n", "\n", "The implementation of the Chi-Square test consists of two steps:\n", "\n", "1 ) Generate the contingency table (see https://en.wikipedia.org/wiki/Contingency_table ) which in our case is a table with the values of **Nij** and **E_nij** for each pair of categorical values.\n", "\n", "2 ) Compute Chi-square value as :\n", " $$\\chi^2 = \\sum_{i=1}^{I}\\sum_{j=1}^{J}\\frac{ (N_{ij} - E[n_{ij}] )^2 }{E[n_{ij}]}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To generate the contingency table we first use CROSS JOIN to form a table with all possible pairs between the two categorical features. Only categories with more than 'MinSampleSize' are considered, 5 is typically used for the Chi-Squared test. The following string performs that operation:" ] }, { "cell_type": "code", "execution_count": 9, "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", "
data1data2
0Infiltrating Lobular CarcinomaYES
1Infiltrating Lobular CarcinomaNO
2Other specifyYES
3Other specifyNO
4Infiltrating Ductal CarcinomaYES
5Infiltrating Ductal CarcinomaNO
6Mixed Histology (please specify)YES
7Mixed Histology (please specify)NO
8Mucinous CarcinomaYES
9Mucinous CarcinomaNO
\n", "
" ], "text/plain": [ " data1 data2\n", "0 Infiltrating Lobular Carcinoma YES\n", "1 Infiltrating Lobular Carcinoma NO\n", "2 Other specify YES\n", "3 Other specify NO\n", "4 Infiltrating Ductal Carcinoma YES\n", "5 Infiltrating Ductal Carcinoma NO\n", "6 Mixed Histology (please specify) YES\n", "7 Mixed Histology (please specify) NO\n", "8 Mucinous Carcinoma YES\n", "9 Mucinous Carcinoma NO" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "query_expected= \"\"\"expected_table AS (\n", "SELECT data1, data2\n", "FROM ( \n", " SELECT data1, SUM(Nij) as Ni \n", " FROM summ_table\n", " GROUP BY data1 ) \n", "CROSS JOIN ( \n", " SELECT data2, SUM(Nij) as Nj\n", " FROM summ_table\n", " GROUP BY data2 )\n", " \n", "WHERE Ni > {0} AND Nj > {0}\n", ")\n", "\"\"\".format(str(MinSampleSize) )\n", "\n", "sql = ( 'WITH\\n' + query_table1 + ',' + query_table2 + ',' + query_summarize + ',' + query_expected +\n", "\"\"\"SELECT * FROM expected_table \n", "\"\"\") \n", "\n", "regulome.runQuery ( bqclient, sql, [] , dryRun=False )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the resulting table has $I * J$ rows. Next, the contingency table is generated by using an \"INNER JOIN\" and filling the missing values of **Nij** with zeros (the IF statement in the query below)." ] }, { "cell_type": "code", "execution_count": 11, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
data1data2NijE_nij
0Infiltrating Ductal CarcinomaNO768705.176471
1Infiltrating Lobular CarcinomaNO118182.420168
2Mixed Histology (please specify)NO2627.226891
3Mucinous CarcinomaNO1715.428571
4Other specifyNO4341.747899
5Infiltrating Ductal CarcinomaYES971.823529
6Infiltrating Lobular CarcinomaYES8318.579832
7Mixed Histology (please specify)YES42.773109
8Mucinous CarcinomaYES01.571429
9Other specifyYES34.252101
\n", "
" ], "text/plain": [ " data1 data2 Nij E_nij\n", "0 Infiltrating Ductal Carcinoma NO 768 705.176471\n", "1 Infiltrating Lobular Carcinoma NO 118 182.420168\n", "2 Mixed Histology (please specify) NO 26 27.226891\n", "3 Mucinous Carcinoma NO 17 15.428571\n", "4 Other specify NO 43 41.747899\n", "5 Infiltrating Ductal Carcinoma YES 9 71.823529\n", "6 Infiltrating Lobular Carcinoma YES 83 18.579832\n", "7 Mixed Histology (please specify) YES 4 2.773109\n", "8 Mucinous Carcinoma YES 0 1.571429\n", "9 Other specify YES 3 4.252101" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "query_contingency = \"\"\"contingency_table AS (\n", "SELECT\n", " T1.data1,\n", " T1.data2,\n", " IF( Nij IS NULL, 0, Nij) as Nij,\n", " (SUM(Nij) OVER (PARTITION BY T1.data1))*(SUM(Nij) OVER (PARTITION BY T1.data2))/ SUM(Nij) OVER () AS E_nij\n", " \n", "FROM\n", " expected_table AS T1\n", "LEFT JOIN\n", " summ_table AS T2\n", "ON \n", " T1.data1 = T2.data1 AND T1.data2 = T2.data2\n", ")\n", "\"\"\"\n", "\n", "sql = ( 'WITH\\n' + query_table1 + ',' + query_table2 + ',' + query_summarize + ',' + query_expected + ',' + query_contingency +\n", "\"\"\"SELECT * FROM contingency_table\n", " ORDER BY data2, data1\n", "\"\"\") \n", "\n", "df_contingency = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n", "df_contingency" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this contingency table, we can use python to compute the Chi-Square statistics from the contingency table. This is used to validate our BigQuery implementation:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chi2, p, dof : 309.39167298726557 1.020302231598449e-65 4\n", "Expected nij : \n", "[[705.17647059 182.42016807 27.22689076 15.42857143 41.74789916]\n", " [ 71.82352941 18.57983193 2.77310924 1.57142857 4.25210084]]\n" ] } ], "source": [ "yes_a = df_contingency[ df_contingency['data2'] == 'NO' ]['Nij'].values\n", "no_a = df_contingency[ df_contingency['data2'] == 'YES' ]['Nij'].values \n", "conting_table = [yes_a , no_a]\n", "\n", "chi2, p, dof, expected_nij = stats.chi2_contingency( conting_table ) \n", "print( \"Chi2, p, dof : \", chi2, p , dof)\n", "print( \"Expected nij : \")\n", "print(expected_nij)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following query string computes the Chi-square and the Cramer's V statistics from the contingency table." ] }, { "cell_type": "code", "execution_count": 13, "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", "
IJNChi2V
0521071309.3916730.537477
\n", "
" ], "text/plain": [ " I J N Chi2 V\n", "0 5 2 1071 309.391673 0.537477" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sql = ( 'WITH\\n' + query_table1 + ',' + query_table2 + ',' + query_summarize + ',' + query_expected + ',' + query_contingency +\n", "\"\"\"\n", "SELECT I, J, N, Chi2,\n", " IF(I > J, SQRT( Chi2 /(N*(J-1))),SQRT(Chi2/(N*(I-1)) ) ) as V\n", "FROM (\n", " SELECT\n", " COUNT( DISTINCT data1 ) as I,\n", " COUNT( DISTINCT data2 ) as J,\n", " SUM(Nij) as N,\n", " SUM( (Nij - E_nij)*(Nij - E_nij) / E_nij ) as Chi2 \n", " FROM contingency_table\n", " )\n", "\"\"\") \n", "\n", "df_chi = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n", "df_chi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The computed Chi2 using BigQuery can be compared to that obtained with python. A large value of Chi2 indicates that the null hypothesis (that data1 and data2 has no association) is rather unlikely. The degrees of freedom ($IJ-I-J-1$) and the Chi2 value are need to compute the p value of the null hypothesis. " ] }, { "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 }