{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Regulome Explorer Spearman Correlation for numerical data\n",
"Check out more notebooks at our ['Regulome Explorer Repository'](https://github.com/isb-cgc/Community-Notebooks/tree/master/RegulomeExplorer)!\n",
"\n",
"\n",
"```\n",
"Title: Regulome Explorer Spearman correlation\n",
"Author: Boris Aguilar\n",
"Created: 01-13-2020\n",
"Purpose: Demonstrate how to compute Spearman correlation coefficients from BigQuery tables\n",
"URL: https://github.com/isb-cgc/Community-Notebooks/blob/master/RegulomeExplorer/RE-SpearmanCorrelation.ipynb\n",
"\n",
"```\n",
"***\n",
"In this notebook, we describe the computation of the Spearman correlation coefficient typically used to estimate the significance of associations between two numerical data types; in this case the numerical data are the gene expression of two user defined genes. To describe the implementation, we used the gene expression table of the [pancancer-atlas dataset](https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/PanCancer-Atlas-Mirror.html). More details of the Spearman correlation can be found the following link: https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient\n",
"\n",
"The Spearman correlation is also implemented in the ['Regulome Explorer notebook'](https://github.com/isb-cgc/Community-Notebooks/blob/master/RegulomeExplorer/RegulomeExplorer-notebook.ipynb) which computes statistical associations between other numerical data types available in the pancancer-atlas dataset, such as copy number variation, protein expression, clinical features, etc."
]
},
{
"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": [
"%matplotlib inline\n",
"from google.cloud import bigquery\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": [
"## User defined Parameters\n",
"The parameters for this implementation are the cancer type and the name of the genes for which gene expression data will be obtained. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cancer_type = 'UCEC'\n",
"gene_expre1 = 'IGJ'\n",
"gene_expre2 = 'ADAM6'\n",
"\n",
"\n",
"bqclient = bigquery.Client()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data from BigQuery tables"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Gene expression data from gene1:** The following query string retrieves the gene expression of the user specified gene1 ('gene_expre1') 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."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"bigquery_table1 = 'pancancer-atlas.Filtered.EBpp_AdjustPANCAN_IlluminaHiSeq_RNASeqV2_genExp_filtered'\n",
"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 `{0}` \n",
" WHERE Study = '{1}' AND Symbol ='{2}' AND normalized_count IS NOT NULL \n",
" GROUP BY \n",
" ParticipantBarcode, symbol\n",
" )\n",
")\n",
"\"\"\".format( bigquery_table1, cancer_type, gene_expre1 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Gene expression data from gene2:** The following string query will retrieve the gene expression of the user specified gene2 ('gene_expre2'). This is very similar to the query used for gene1."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"bigquery_table2 = 'pancancer-atlas.Filtered.EBpp_AdjustPANCAN_IlluminaHiSeq_RNASeqV2_genExp_filtered'\n",
"query_table2 = \"\"\"table2 AS (\n",
"SELECT Symbol, data, ParticipantBarcode\n",
"FROM ( \n",
" SELECT \n",
" Symbol AS symbol, AVG( LOG10( normalized_count + 1 )) AS data, ParticipantBarcode\n",
" FROM `{0}` \n",
" WHERE Study = '{1}' AND Symbol ='{2}' AND normalized_count IS NOT NULL \n",
" GROUP BY \n",
" ParticipantBarcode, symbol\n",
" )\n",
")\n",
"\"\"\".format( bigquery_table2, cancer_type, gene_expre2 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point we can take a look at the combined data (Gene expression of the two genes) by using a simple INNER JOIN command. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" in runQuery ... \n",
" this query processed 7757877633 bytes \n",
" Approx. elpased time : 2798 miliseconds \n"
]
},
{
"data": {
"text/html": [
"
\n",
"
\n",
" \n",
" \n",
" | \n",
" data1 | \n",
" data2 | \n",
" ParticipantBarcode | \n",
"
\n",
" \n",
" \n",
" \n",
" 1 | \n",
" 2.819599 | \n",
" 3.831652 | \n",
" TCGA-AX-A2HD | \n",
"
\n",
" \n",
" 2 | \n",
" 2.180091 | \n",
" 3.760516 | \n",
" TCGA-EO-A1Y7 | \n",
"
\n",
" \n",
" 3 | \n",
" 3.027851 | \n",
" 4.325743 | \n",
" TCGA-DF-A2KR | \n",
"
\n",
" \n",
" 4 | \n",
" 3.159486 | \n",
" 4.401590 | \n",
" TCGA-QS-A5YQ | \n",
"
\n",
" \n",
" 5 | \n",
" 2.006008 | \n",
" 3.369119 | \n",
" TCGA-QF-A5YT | \n",
"
\n",
" \n",
" 6 | \n",
" 3.278023 | \n",
" 4.434266 | \n",
" TCGA-AX-A1CF | \n",
"
\n",
" \n",
" 7 | \n",
" 3.148507 | \n",
" 4.542022 | \n",
" TCGA-A5-A2K4 | \n",
"
\n",
" \n",
" 8 | \n",
" 2.778011 | \n",
" 3.745502 | \n",
" TCGA-BG-A0W2 | \n",
"
\n",
" \n",
" 9 | \n",
" 3.720214 | \n",
" 4.769811 | \n",
" TCGA-AJ-A2QK | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" data1 data2 ParticipantBarcode\n",
"1 2.819599 3.831652 TCGA-AX-A2HD\n",
"2 2.180091 3.760516 TCGA-EO-A1Y7\n",
"3 3.027851 4.325743 TCGA-DF-A2KR\n",
"4 3.159486 4.401590 TCGA-QS-A5YQ\n",
"5 2.006008 3.369119 TCGA-QF-A5YT\n",
"6 3.278023 4.434266 TCGA-AX-A1CF\n",
"7 3.148507 4.542022 TCGA-A5-A2K4\n",
"8 2.778011 3.745502 TCGA-BG-A0W2\n",
"9 3.720214 4.769811 TCGA-AJ-A2QK"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sql_data = 'WITH\\n'+ query_table1 +','+ query_table2 \n",
"\n",
"sql = (sql_data + \n",
"\"\"\"\n",
"SELECT \n",
" n1.data as data1, \n",
" n2.data as data2, \n",
" n1.ParticipantBarcode\n",
"FROM\n",
" table1 n1 \n",
"INNER 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 of the two user defined genes, we can use a scatter plot:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX10VOd957/PjEYjjWQQkS3zJmHAxlg22PgF8EKySmTc\nNLGTutvSZNvjdndPa067J628p7TbPT2Gs93u1t612u321OxJm9Q9aZs0W5I6SesQYiWFNAYb26HG\nGAw2SAIboyCM0OuMnv3jNw/Pnas7M3c0c+fOy/dzDmc0M/cV7K9+9/v8XpTWGoQQQspPJOwLIISQ\neoUCTAghIUEBJoSQkKAAE0JISFCACSEkJCjAhBASEg1Bn0Ap9Q6AywBmAcxorTcGfU5CCKkGAhdg\niPD2aK0vleFchBBSNZTDglBlOg8hhFQV5RBGDWCfUuqwUuqXy3A+QgipCsphQWzRWp9XSt0AEeI3\ntNYHnBsopfQTTzxx7X1PTw96enrKcGmlZ2BgoGqv3Ytaup9auheA91NJDAwMYGBg4Nr73bt392mt\n/zDffoELsNb6fPr1faXUXgAbARxwb7dr166gL6UsVPN/RF7U0v3U0r0AvJ9Kwh007t69u83PfoFa\nEEqphFKqNf1zC4AHAfxLkOckhJBqIegI+EYAe5VSOn2uL2mtvx3wOQkhpCoIVIC11m8DuCvIc1Qa\n1foIlY1aup9auheA91PhDPjZSFVCP2CllK6E6yCEkBKh/GzE/FxCCAkJCjAhhIQEBZgQQkKCAkwI\nISFBASaEkJCgABNCSEhQgAkhJCQowIQQEhIUYEIICQkKMCF1zvAwsH070Nsrr8PDYV9R/cBSZELq\nnO3bgVOnAKUArYGbbwa+/OWwr6rqYSkyISQ/IyMivoC8XrwY7vXUExRgQuqc9naJfAF5bW8P93rq\nCQowIXVOf7/YDgsXAqtXy3tSHugBE0JI6aEHTAghlQwFmBBCQoICTAghIUEBJoSQkKAAE0JISFCA\nCSEkJCjAhBASEhRgQggJCQowIYSEBAWYEEJCggJMCCEhQQEmhJCQoAATUsHU6rSKWr2vQmE3NEIq\nmFqdVlGr9+WA3dAIqXZqdVpFrd5XoVCACalganVaRa3eV6FQgAmpYMKeVhGUVxv2fVUK9IAJIdcY\nHgb6+sQiaG8HJiaAc+fC8Wrd19LfDyxbVp5zlwBfHjAFmBByDffi2OAg0NVlv1+4ENi/P5xrqbKF\nOi7CEUIKw704BoTn1dbDQl1D2BdACKkc2tuB0VEbdW7cCLS0iPgZGyCsa6nFhTpaEISQawwPA48/\nnim4YfmulXQt84AeMCGEhAQ9YEIIqWQowIQQEhIUYEIICQkKMCEkA3YqKx9chCOkDslVZVblBRCV\nArMgCCHe5BLZ3l7g8mX5eXoaePddYN26qkwFCxNmQRBSD3hZBvlshFxVZs5OZYODwMyMCPKpU5KX\nS0oHBZiQKqevT8TRKZJenznJ1Q7S2aksFrO9IGq1HDhMKMCEVDle0Wy+Pgq52kEuWyZ2xP79QE8P\n0JBuWFCr5cBhwl4QhFQ52Xom5OqjYEQ2H/39c8uBSengIhwhVY5XzwSg+D4KVd6PN2yYBUEImT9M\nRysKZkEQUkuUu0CiHvrxhk1ZBFgpFVFKHVFK/X05zkdILZIvs2G+ZBN2Ds4MnnItwv06gGMAFpTp\nfITUHEFFpEbYk0ngyBFgwwZg0yZgfFzygAF5zwW40hN4BKyUWg7gEwA+H/S5CKllskWkxVoTRtgH\nB6XybXISOHBAxLirC+jsBBIJLsAFQTksiH4AvwmAq2yEuChEPN25uzt3yj4bNgD79klEPB9rwgh7\nMimv0SiQSskfgP5vkARqQSilPgngPa31q0qpHuRYGdy1a9e1n3t6etDT0xPkpRFSEZjHf6Ukb/fx\nx7NnGrhzd02WwuSkiOXgILBqVeFiaXJ9TdlxZ6e1HgD6v34YGBjAwMDAtfe7d+/u0VoPZN0hTaBp\naEqp3wfwCwCSAJoBXAfg77TWj7q2YxoaqUucjW+Awsa+m31PnQKmpqRi7ZZbJDr+ylcKvxZnPnEi\nIZ+NjzMHeJ5UVh6wUupfA/hPWutPeXxHASZ1iTvXthDxNPvOzEjEGovJFGOlKJwVAPOACal0cvVk\n8LvvDTcA27YBr7wikeu5c+xeVi2wEo6QGqIYS4OUFEbAhNQbLJ6oLijAhNQQ87E0OAMuPGhBEFLn\nsOlOINCCIKSemG8ky6Y74UEBJqQCKIUNMN9mPfSNw4MCTEgFUIh4ZhPr+UayxaTCkeLgSCJCKoBC\nxDNb+XK20UT58DueiJQeRsCEVABuG6C5ObslkU2sGclWH8yCIKQCcM91m5iQirZkEjh7VsqMe3pE\nVJ0RcKHly0FdO2fHzaGyekHkvAgKMKlh5iNQW7cCx44BV66IyDY1AbfeKhHu008XP3CzlPhNY6sz\noaYAE1IJ5BMoL2HavFneT00Bs7MSAa9fX5mlxX7Ln+ss35h5wIRUAvkW2LwyILq6JOqNRKRBeiyW\nfWEt7Eo2v2lszDeeCwWYkIDJJ1BOYZqZAV54ATh5Ura96Sb5bnoaeOMN4O2354psUMM6/eJ38Y/5\nxnOhBUFIwLgX2HbuBJ580loO4+PA+fMitKdOyT5mKsX4uLSY1FpEuKkJWLky8/HdWACmLzAAPPDA\n3POE7bm6/x7Cvp6AoQdMSCXi9kI/9CHg+HHg6lX5s2oV0NIi254+Le9PnJCxQ9EosGZNps9qjvf2\n2zKeqLHRNmVPJMTOaGioec+10qAHTEgxlMpb9TqO0wt97TWgo0OENpEA3n1XvtNahFhrEVAzMNP9\n+G4sAACIx220nEzK6+AgPddKhREwIVko1ar99u0S4Q4NiU2QSsmxGhvluIODEqUC8v3588C6dSLG\n4+PA0aOSCdHYKPbDsmXej+/mek+elONMTcnn0Shw223A2rXh5gvXGYyACSmGUq3aj4yI+JqUskgE\nuHTJLlpt3GgXpxoapOBi/36phhsdFXFesQK45x4R3ZERWXhzR+QmEm5qkui3oUHOBcj5WBlXeVCA\nCclCIav2ueyK9naJSM1xYjHxcffvl4h0zx6bRbB0qUS9vb3AwID8bKLab39bMiGyZTuYng6vvAIs\nWiRCvGAB0N0t53NGzMPDwMMPi/XR0QE89BAbsYcBBZiQLJiIMh4HLlyQ0uBsXvCOHcC+fcDhw/K6\nY0fmcdraxAqIx4Hly+eKuRH6V18FzpyxWQ1vvSWRcyolf8y5lcou+suWSRR9yy3iKzc0zD1fXx9w\n4IAs+o2NAQcPcoBnGNADJiQPfrzgjg4RM0NLi4i24aWXgEcekW1aWiTq/eIXxU44cUIEOh4H3nxT\nbINIRATYZDI0Nop9MTsr0azW4hVPTNjsiK1bgeeek/PlS/nq7ZVfFqmUvI9Ggfvuq7wquyrGlwfM\ndpSE5MGvF6y1FWk3Tz4pIm2+f+wx+3501KafxWLS/yEWk+9MFdyqVZLRcOmSCPXZs8D779vFuWQS\nePFFe758LSbb2+XYqZTNrshWZVdH/RvKDi0IQvLgxwveuFGE0dgMGzdmfu8W8atX7ftYTAQUEHsi\nEhHbIB6XaDsWE3947Vrghz8UAezosNczM2OP5Zf+fomYW1qA1lZgyxbvRbqwq+xqHUbAhOShv3/u\n47zBRIiXLomQLV4sebyjo2JdmIjR3Szd5PcqJaJrsiLa2yXrwVTGaS0pac70MSPmLS3i32otUfCm\nTf7vadkya1c4cUe8586xf0OQ0AMmJAf5HsHd/vCFC5lWg/GLvcqRn3rK26P18m8Bex3GM45E5vYK\nLtYeyHc/YfceriJYikxIseRbgHO3YjSlw4ZStY90XofxgtesKb0v674fk7VRJ/0bSgkX4QgplnwL\ncLmsBbdfnCuazhdpO6+jsdHmEZca9/1wXlywcBGOkBzkW4DbuVMe00+fltc9e6SY4uxZKTGemLD5\nubkWtB57LHsesZ/rKIRcRSOcK1deaEGQuiZf5Jkvn9bLotDa27bINTmio0MW1Mw+ra2ZecSlbOVY\nZ5MpwoIeMCH5KFaMvEQVsJ85m+ucOCElwqYJj3NByxRyzM7KPpEI8FM/FYznmu0XAXN+Swqb8RCS\nj2Ib7nhZA83NIuonTshgzclJEby2Nptu5uz5sH27zHszBRVaS/VbUHm32ewM5vyWHwowqWuK9Va9\nPFMj6KmUCKppqKOUXTxrbhaf+JVXgG9+U/oymF6/kYjsd/q0dFErNdl8Xs5sKz+0IEhdU+oxOcPD\nwIYNEvWaPg1KSTpXPC4VZ4kE8J3viOcbjYrlMDsrEfKVK7JPJCKfNTVJJ7RyWAFuO4Y5v0VBD5iQ\nQjE+6Llz0pWsqyt7A3S3Z7pzpzTcOX9evp+dtf0cmprkz6ZNcuzTp+VR3whtJCICbTqfAbJvIgF8\n8pPlWSSrs5ltQUMBJqRQ3PPVzBDMpUvFNhgZsUMyDx2S6NXMXLtwQTzeiQkRVUCE9Y475PvVq2X/\ny5elmOL4cRvtRqOyTyplhVsp6ed7993sUlaFsBCDkEIxPmgymfn64osyqVgp4MgR2db06D1zRkTU\n2AfGx52dlaj2+uttRNnXJ0I9NGRbTHZ3A++9B3zwgR2maYTZq3cwqR0owIQ4MJVgDQ0ivmZhDMhc\nXAPku2RS0sdM+0jng1wsJn6wM3rt7wc2b7bHjkZlMa6nRyLnc+fkO9Pj4bbbWAxRyzALghAH/f1i\nN0QiNgpdskS8WyOuRjg7O0UkUymxIuJxeR+J2PaRbg912TLJhFi5UiLdK1fEtjh+XL6/+WaJmLdt\nkwyJr3yldD5sqaY8k9JBD5gQF17FGU8/bReoEgnZbnxccn0/+MD25G1oEN92xYrsi3jbt0u58eXL\nmV6vO1ouBD9FFKyAKyssxCBkPmTLhx0fB157TYZj/uAHEumOj4tlYAQ4Hpem6W1tkmZ29Cjw/PPA\nvffayPOXfknsBq3FA47FZH+311tIxOqniIJ5vpUHI2BCXHjlw46Pi5CayRUmvQywpcXOPN+vfU2s\nicZGyXhQCrjrLttjd2xMRNgs8i1ebKddZLuOXBFrrj4Tue6Leb6BwQiYkPngVSl26JBdfANEwIzo\nNjSIyI6OAt/6FvB3f2eLK8w4etPH9/RpmZhhiiwSCekP4RZfoLCI1U9FHzudVR7MgiB1Rza/1P35\ns89miqJZmDMoJfYBkCl+bqGORERoBwel0EIpEeamJhHC1au9F9rcvXlzpaPlGptkYG/fyoMWBKk7\nsj3a53rkf/hh4Pvfl5QzYy185CMiogcPZuYAA3a0/KJFdlbbwYPyeuONkvcLAA88kL3ijJVpVQ0L\nMQjxItujfa5H/meemTvT7cknZZ+WFhFfMyATEIshGhXxfeYZ+WzzZolo33tPCizWrs3twTJirX0Y\nAZO6wx3pmjLjgYHM0uJci1TOY0xNiRCbHhBO2tpk/HtzM/DGGxLVJpPiw3r5vqRm4CIcIV64F6PM\nBIvFi+X74WHJVDBpYF7pX85oOR6XirVFi2yOsGF6WsqYz52z4tvQIHnCXs19WChRXzACJnWLWXT7\nznfEMjCZDdPTUq2mVPax7yYCnpmRxbVYTCyI2VmJiA2m0U5Dg7UlzILa4GDm9bBQoqZgBExqC78R\not/tTPECIItoV69KhJpKifC+8YZ4thcvSnrZ5s32WCaKfvddeb9kiUS1RnCVkp8BsR+UsgIcj4vN\n4YaFEvUHBZhUDX5H5vjdzgheZ6eIpdY2f/fyZZtOprW0phwdtccyC2Tr1omNEYvJYtz990sJ8s/+\nrFgcCxcCN90kotvUJJH1qlXZ085KNfmYVAcUYFI1+I0Q/W5nBC8WkynEWksEbHr5mldAxDganXus\nbKKptdgO5rvly60gZyuCYKFE/UEPmFQNfktpvbYzvXidxRcAsGOHLJLNzsp3uWhtlayGNWsyj+HO\n1TURuNMfdnvIpOYJfyKGUioO4PsAGiE5x1/VWu/22I4CTPKSrTDBLKYND4t3u3ixeLMrVkiKmVMU\ncxVfvPzy3HOaSRWmCi4el/eNjbaEGMgU9+HhzIU4r74MpOYJvxBDaz2llPqo1npcKRUFcFAp9Q9a\n60NBnpfUJtkKE4y4nj4twnf1qvTbdW6fr/jC9GxwY8TXjAkyFsXsrPWETRqbUvLZhQsizkNDcty2\nNhHlfNGvn5aSpLYI3APWWo+nf4xDBJ+hLikpRkTNBGLTYczp17q92uZmiX6PHrXi7aax0S7OATaT\nwUTRsZicY2TEjpF/800R4OPHJbNCKRHgbAuBTvwuHpLaIfBSZKVUBMDLAFYD+BOt9eGgz0nqC9O0\nxowIMotfiYRNQ3v7bcnvjUSAjRtFGE+dkvSxs2clao5ERHSNiCcStnnO1JRkOZgWkg0NmfPajhyR\n409OyvtIRAS6oUFsCz8pZUxDqz8CF2Ct9SyADUqpBQC+ppTq1lofc2+3a9euaz/39PSgp6cn6Esj\nVYzzcT2REK+3qUnE1Hi/4+OSy2umDzc0ALfeKkJ68WKm2JlXE90mkzLpQmv57PrrRbh//OPMc5iF\nuA0bxG4wxRwmCyKZ9E4p87IbCul+RiqLgYEBDAwMXHu/e/fuHq31QNYd0pQ1C0Ip9bsArmqtn3Z9\nzkU4UhB+qsZ6eyUyNaN/jCi2tooIG592ctLm/xq/15mCFokAH/uYbJ/NnzVjhpyibbqheWVA5Bt7\nRA+46gl/EU4pdT2AGa31ZaVUM4BtAP5HkOck9YGfx/X2djtY0wiqiUjb2mQYZiplLYVUSrxhI8hG\nHCMRGUV0440i2KkU8M//nNlMZ+dO4Otft+eOx0XoX3nFW0S9rp/dz+qPoBfhlgB4QSn1KoAXATyv\ntf5WwOckdYBZVJuelkjy6NG5Zcf9/ZICZqZWAOITGyG9elWE1vTubWgQO8MIrzMKBmSfqSkRYGdV\nHCCtKRMJsUGamuScPT3ZI1hWvREgYAHWWh/VWt+ttb5La71ea/3fgjwfqQ+Gh8XfHRyUrIPZWfFo\n9+0TL9YI8bJlEqU+9BDw4Q/L+7VrpR+vydN192iYnJw70eKGG8T/NVkWgM2AMNdjWlmayDkWy13J\nxqo3ArASjlQhTv/0xAkb3RpRNUMwb7wR2LtXJhIDtpBj3z5539kpIj49LYKZSklUHIvZfN9YzKao\nmYbqsVhmQ3Xj/05P2zlx27Zx4GWdE34lnF8owCQX7owBZ6XZqVM27zeZtGlgZkS8V9vHhx8GDhyw\nEe3UlIhqLCaLaIDsC8hi3YUL9jq8Fsl6e+WzwUG5hqam7N4vqRvCX4QjpBiM8JrH+xtvlKyGiQlZ\nLOvslEh0dFREdGZGvjPpZEpJROtmfFx695oR85GIpKfFYsDJk7KPaSnZ0CACa8YQecUJJn1s1Srb\ne6KU4htkhRyr78KFETCpWIzVcPKkiKVJFQNENBculPQwZ1Ocr31NtjURrFcE3NEhgmPydQER2u5u\nW2K8YoWIPSAR7fLlkjXR0TE39S3o4ZlBNmpnE/jAYARMqhuTqmVSxIzVYCrWIhHxWZ3NeBYsEBsh\nmZSFub1750Z5RsSdv/O1lgW9BQtspVwkIpbG1JREzM3NEoUDmalvQaePBVkhx+q7cKEAk8CZ72Nu\nc7NEoTMz1ts1UaszU8HdjKelRZrx3HyzLMA5o7zRUSvezmjapKHNzAAvvCBWhlOgZ2fluxMnpMF6\nLFa+1LEgK+RYfRcubMhOAsOMBtqwQbIE3n+/sCYzJjIzAmlG/MzOSsbB+vU2BezkSTsW3tgVL7wg\n/q3xkM0xV64EPv5xsSmiUSmYiEbF+x0dFY84myM2MSGtLsuZOhZkyhrT4cKFHjAJDLeH29QkC1V+\n++P29opoDw6KreAujFi8WCLc739fxNNExQsWWLFevdrOfTMTkE0jd6d3OzAw9/heXHcdcN997O9L\n8sKhnCRcnB6uSRMr5DHXLKCZjmQGk+Fw8SLw/PMSsZo+DoBNKTODL7u65L2J8nbulF8Ojz4q1/Ps\nsxIFqzz/y5iBm7mun6PlSSEwAiaBYSLgZDL7ePdcDA+LfTE5KeJnphY7cU6saGwUW+Huu0Uknav7\nS5eKpzwyIj5uW5tsa1b+f/ADibbNVGSTZXHlipwnEpHod/NmYM+e7NfPrAKShoUYJFxKkZ7lFLST\nJ22hhMFE186mOw0NYkOYRjrr1wPHjolHHI1KRG3sEECE9g/+AHjkERH5lpbMCrpC6O2V7msGjiOq\nW5iGRsLFT3pWvgyJ/n4r4vG4RLHT09bvNRkNY2PWw52Zkb69CxdKPu/rr0t068ycmJoSYTdN1e+9\nd26+8HxgVgEpBEbAJFSy9cX1EmVnz4VUypYSe81zi0Qkkl2zBvjRj+Zuo5R8f911YllMTOSP0v2k\n0wVdlEGqBloQJDz85v46Mx1SKbEGNm4Ezp+fO9Z9505g924ZIw/Idt/+dnYBXrBAbIbXX8+sogNE\nfNeuFW/a2ed34cLMPr9O6O+SAmAWBAkPvwMmnZkOqZQI5aFDInLm88lJOcZTTwHPPCNifMst0iTd\nS3wbG6VkuLtbsh4+/OHMPGKzYGd+5w8NSQQ8Pg6cOycLbV7ZC6waI6WGAkwCwS1WQ0Pe6Vn9/ZnD\nK03qmKl201oE+ORJKax47DER42PH5BxenD4tEfTBg3L8eFwi21hMRgR1dIgwX7ggUfHYmNgaZuKx\nu9m6gU3USamhAJNAcIvV2bNzI2JjUwCSndDRIdvNzoo4NjZKhJtKiUheuiSjgJwj6COO/4KVArZs\nybQP+vpEjG+6Cbj9dunTe/68FHGMjYm4m/xkExk7m607YdUYKTX0gEkguBejzp2z/RwAETGTq2vy\nhMfHZazPkiUikpGIiC5gezWYcfMm68GknkUiUkyxdSvw3HP2PNnSwjo6bKtKE2W3tMxtth7U3w1b\nQNY8TEMj4eFOQXMvYLW3W5siFpOI8vRpWTQ7fVosgWhU/pgpE4AtxIjFbPRq9m9pERF3kistzEwv\nBoAPfcg2VjeiWCh+hdX44067g4t59QkjYFIWvNKznEKktUS9ExNSfaa1COrEhJ1anEiI2N50k+18\nZsYJNTVZC6OpyRZTLFli7Y4zZ8RjXrZMRPLIEbEyolGxLr7xjezX7kdY/WZJsFijLmAaGqls3KJs\nRNHYC6aPxOysiOrChcCdd4pQmwY/sZgIKCARsRmwaaJdU1zhFselS0XQ/eTrllpY3cczzYFITUEL\nglQ2bpuit1fEyOT/XrkiItnZaZvp7Nkjoj04KNuZ71avBr7zncxpx1evAi+9JCXG774rnzmtCqdX\nnAu/6Wd+q+Cc1X3ztTtIbZAzAlZK/TGAbBtMATgF4Eta6ytFXQQj4JrHz2O8OzK8cMGOAJqaEnFb\ns8bOZ3vqqUwR27zZiqVTAEdGbPVcLAasW1dY1Ok3Yi20Co6LcTVN8RaEUuoXc+zbAOB2AOu01tsK\nu7Y556EAVxF+S3Kd20xMSCZErsd4t4AZkR0akinDkYjNUujunru/iXadDXUefFDEW2s7rv6RR2zU\n6ec+duzIrL7L1Q2tEFhZV9MUb0Forf8i71mU+pbfKyLhUOpIy88qvnubwUFbZOF+jHfOdDt7Vhro\nALKA9vTTEtmansDJpGzX0TH3ukxDHXO83/ot6wtHIiJyxi8u5D7OnZNrNwuDpYpSWVlHcgqwUmor\ngFVa62fT778K4EPpr39Pa/1drfUnAr5GUiSlTnvyIxzubQA7Lsj0Bt6+PTMbwmQ2XL0qEfPjj8tx\nhobscVMpKZZobpb9vX6pOO+3q0vOl0yKCK9aZQtB5nMfpRRJdk4j+SrhdgN4yfH+VgC/CWAXgJ0B\nXRMpMcWIiNeEBz8lue5tNm6UR+zz5+WzxYvnCqGpbjPTj4eHpfzYzcKF8n22XhPO+21tBe6/X/pB\n3HmnRLDm78DPfSQScvwTJ+Q1kfD/d5cPVtaRfAK8QGt9zPH+pNb6Za319wFcF+B1kRJSjIh4NdXx\nIxzubfbskajbLIA1NmYK4dSURL0TE3YU/Jkzc+e0KSXdysbHs/9S8RJWr8/83EeQSxMmC2T/flnU\n4wJc/ZEvDa3N+UZr/dOOtzeW/nJIEBQjIl7Rs59G6+5tTCR99KhdFAMkQu3ull8OZppFLCaP5l1d\nkj7mHEO0eLEcO9fju0nzMsUXTU3ShGfpUhFup2XhtRDoXHRLpSTVrbFR3rsr7QgphnwCfFwp9Umt\n9TedHyqlHgLwZnCXRUrJxIREeIZCRKRUPqWJpJcskU5mZkHr/feB731PtjHerpmc3N4u4vfOO7ZZ\nzuc/L9vu3JmZ8bBnjz2XEdbt2+Xep6bE+rj55vy5v319wIEDkjGhte0jbCYqF+PTMu2MuMknwH0A\nvqmU+hkAR9Kf3QPgXwF4KMgLI6WjGBEtVdGAs+9Dc7OI2vh4ZnQ7OyuC9+qrkq3wkY+IwEYism8i\nAfzpnwKf+ATw5JM2R1hrSVdzT9IYHi7c+x4ZsaOLlJIeFKYIpNiiCfaAIG7ypaG9pZRaD+DnITm/\nAPB9ADu01pPZ9ySVRDEi6sdu8IPzl0A0KpGpc0YbYP1e0zz9pZdkO2cjHmMNeFkjboEzHdic/R78\nXGc0ansRNzZKA/hSlAoz7Yy4YS8IUjDzeZR2FlkkEjIG/oMP5o6Zj0TEswWsUDsX4hYuBB54ABgY\nkGi5q0t849Wr5XouXxb7YHBQekoAIqKNjXNbVWa7TqcHvG6dXK/bO54P7AFRV5SkEu4KvEuRFQCt\ntV4wv2ubcx4KcBWRr4LLb9nx8ePAv/zL3EXCeNymo5lIFLCVcLffPndenFc+8cyM7ZS2fLks6K1b\nZ6vsnnwy8xqBudft7thWTLUaB3bWFeyGVq8EvdizdasspJk2kd3dsnBlMALtLLowImmuw4jR3r12\nrpsROTMJQyn7WSQiE4xTKZkHZ4jHbXvJSEQE/eJF2a+pydoPpgLOLKY5+0wYYdV6rthevMjWkWRe\nsBtavRL0Ys+ZM+Ktmij1zJnM743X+c47soimNfD881Iq3N1tI1CtRVR//GO70NbaKgJsmq1PTMgx\nIxER2tEOz0bwAAAdLElEQVRR20h9akp+EfzoRyLyqZSdjDE1ZY85OWknaZio+OpVbz/W/Rmr1UiQ\ncCZcDryqwKqBoBd7urokuoxGJQKemPCulBsft9VtY2M2mjx1SlLITp0CbrhBvjf9GpYskXNEo7b/\nA2Cj08WLJXqNx0UYIxE5xwcfSPvKsTHZpqnJjqZfs0ZS1Yw3bFLgSlWsQch8YQScg2pNGwo6alu2\nzEbAp06JAF6+bP+OTNbF0aMipI2Ntk8vYHv1KgW89558n0qJaI6NSdnymTPiEZvijOZmuRetxTpY\nvly2P3JExBeQ72Zn5ZgrV2YuzHV2yjHHxkSE77xT9ncurgHeHm01/JuT6oQRcA6qNW0o6KjNHD8e\nl+h3ZkYWvpJJWyn39NPA9dfL9sbPBWw5dDxuR89Ho2JFrFkjf371V62way3d0cykCwDXxtwfOybi\nOztrMyUSCdnvwgWJxk+ckMi3sVEi4pYWieB//GPZ1l0G7LUUUa1PQqTyoQDnwE+zlkok6B4D5viL\nFonYjY+LEL7zjv076uuTn1tb5b0RSSOUd9whIt7UJOLY2Wn/jh97TLZpahKxPXsWaGuTqBewY+7H\nxzNT1IznOzMj0fgbb4jgXrokv4xisewtMc01ezX4yfY5IcVCAc4B/b/cHDokUakZDT85af+ORkYk\nyo1EbB+F2VmxHsbGRBS//GVptP7ggxItm7/jDz6QqNV4wE1N0oCnu9v+W6xYYbMjTHQN2MW32Vn5\n+cIFiar375dMjIa06eb1CzXbE4/z85kZ6dCWKxpmxEz8Qg84B/T/8uNM8XI2Kzc+dDIpYuh8tE8m\nxT4AvP+OZ2ZsRkMqJe+9xty//rp9bxbcYjHZ3oj3zIwV2nwVgdm8c+fnZsin0/PO14y+WtYOSPlh\nBEzmzfr1EqlOTsrr+vX2O/P00NQ0t9oNEMHMxu23i5ACtvDCTX+/5CMvWiQR9oIFtsfE7KydltzW\nZoU2nzWT7YnH+Xk+GwOo3rUDUn4YAZN509wsHq9ZSGtutt8ZsRseBjZskK5nTtx9fp2sXGnH0mst\n752YQpPxceBjH8usWMtV/JGPbE88zs9NkQmQuxk9c4eJH1gJR+ZNb2/2KjFnNd6JExIFmj7ASgH3\n3AMcPux93Hwlu16l0OWqWPNTTsySYwJWwhE/FFO2nCvSMxHpzIxdVDNVal1d2aNaP9cxMmJ7QaRS\n8rppU3mizvk0oyckG/SA65xiUqxyZYkYH3RwUMSyqUlyfRsbgbVr5y6AFXId7e1y3Kkpu0gHMGOF\nVB+MgOucYhaMckV6iYStUjM9HlavFoH0asHo9zqGh6X4w/SYSCQkoh4fz99qkpBKgwJc58x3wcht\nGbjbO5qxRyYnF8g8vnv/RMLfdfT1SaP11lbJvjD9KLjQRaoRLsLVOfNdMHK3nBwft9FoQ4NYBF1d\n4v0ODck+Dzxgj+9cSJuakvObCRmbNgHPPON9HWbhz3jA7uOGDee+kTRchCNz8RIIr+nF+QTE6fFO\nT4sQm25jq1bJNqa3r2mM47QenAtpY2OybXe3pI8lEtlFy0TssZj3ccOGRRikELgIV2fkW+zyuxhm\n+mQkk7Yk2Pl+0yZbiGEa4zjLcp0LaaZSbmgovw9dyvLwIEqGWYRBCiFQAVZKLVdKfVcp9bpS6qhS\n6nNBno/kJ59A+BUQI4RmekVDg7w2NoowPvOMRH5Ll0r7yKmpTEHv75cotqFBfNx43Ip3Lj+3lI2G\ngmiyU60NnEg4BB0BJwE8rrW+HcD9AH5NKbU24HOSHOQTCL8CYoTw/vtlQSwWk9f7788URregDw1J\ntPnooyK6K1YAt90mPzc1lTeFLIholQ2cSCEE6gFrrd8F8G765zGl1BsAlgE4HuR5SXbyNaTp78+c\nCrx0qTyaZ4s0x8dFaIyfe/CgCKzJijh61DZjn52VhbYrVyRSXrRIuqKtWQPcemv5F6yCKBlmEQYp\nhLJlQSilbgIwAOAOrfWY6ztmQVQIw8PA5s3yWB6NSg/e7u7somKyGd56SxbTTLVbU5OIdzIpXctM\nQ3Wt5fubb5b93SXD2bIIgsguYMkwCZDKmYqslGqFiO9/1Vp/3eN7/cQTT1x739PTg56ensCvi8xl\n+3bgm9+0i2MNDTIiyDn12IlzunEqZSddTE2J0DY02KkVzc12lNFdd8l27iyGbCPvs33uvA6mf5Gw\nGBgYwMDAwLX3u3fv/qjWeiDrDmkCF2ClVAOAbwD4B631H2XZhhFwGcklVr29wD/9k0SuRuyWLbM5\nt9n237BBqtMAWyLc3CzHGB+3C22mXeTWrd5Cma3BT67GP0B+gSakzFRMHvCfAziWTXxJ+cmVq9re\nLgtqRiijUdv/1vDYY+L1mjaUO3ZI2tmBA5ILbJqpJ5NyrMZGEV9Aft66NXvZsJ+m6IVMsyCkkgk6\nDW0LgJ8H8DGl1CtKqSNKqY8HeU6Sn1xi1d8PfOhDUgxhmue89JKkkj38sES/hw7ZKHdqShbsnnlG\nRgu1tIjt0NhoXx98EPjJnwTuuw/Ytk22zYZXFoHp/zA4KFV3S5Z4T7Ng+hepNliKXGP48ULdj+tu\nH9b4ui+8IFkKDQ12ttvWrcDzz9sOZJGIZDMYEe/tlZ8HByUCbmqSuW/F+LF+7AUuqJEKo3IW4fJe\nBAW4ZJRKrMwkCxMtm4i2qUlE2TlmaMkSaZDjdf5SlArn838JqUB8CTBLkWuMfF7o8LB4ti+8IDm6\npmuZm74+G+XOzoq3a4ZvNjhWDpSSdDNDEIUItBdIrUIBrjHyiVVfnyyWmfHwBw96l+COjACdnZJK\nZkbPb90qKWkm+jVR7tGjEvm+9JIc368N4LcXA6vLSK1CC6LGyGcv9PbKLDbT+jEalcUx9yP9Qw9l\nZjps2QJ84xty/DvvFEsglRIRbm6WSrYLF2Sxzsv+8PKmndkYTB0jNUbFpKGRMpKvFLa9XQQ1lRLR\ni0ZFADs7JSpuaZGiCuX6z8f5vrlZsh8mJ8WOiMXk+6tXs9sfXqlvTB0j9Q4tiDIRROvD+dDfL1ZC\nS4vYC1u2AMeOiRhOTcnrI4/YHg9r1sir8Yr7+oC2NhFdQOyIzk4R83hcRPbECXl1jqn3Elt6u6Te\noQCXiSBaH86HZcukCOLCBUkPSySA99+XRTYzMv7qVSuOMzNyva+9JkL7/PPS0cz0iFi0CLjhBhHp\nO+7IPJczavYSW+PtZusZTEitQwEuE4U+bpcjYt6xA9i3T+wIU1ShtUTHRhzPn5dtp6flHiYmZLuh\nIbEfenpsb97ZWe+oGchcSFu6VL579FE5X1ubd89gQmodCnCZKPRxuxwR84svirA2NmZe19691kte\nt07E1DA7K97v6KiI886d9rtc9+hspN7cLPuaezt0iF4wqU8owGWi0FSqci1QmYW45mbg+uulgu3e\ne+33RlSjUVt+bEYQTUwATz2VaR1cuJC/sbr73sx1mFd6waReYBZEmSi0UXcQzcLdbNyYmWq2cePc\nbUwD96YmyfM1XdJiMdnv4sXMDIeODn+ZGM57W78eePNNm4WxZ0/p75WQSoR5wBVKUL0NnPm4znaR\nJrXstddku40bRQid59y+XTzjqSlbnrxtmxzz2DEr5N3dIux+721iQkqZmQ9Magj2gqh1Cm1Cnmva\nhRFXkw0Rj0sXM3fTc+e4ok2bpLPZ5s3WVjDRurN/cD7Y64HUIOwFUesUulDX1yeP/ibjYXjYessj\nI7Y6TilrLzhxprBduCDi29cnEy9MeXIyKe8LydxgPjCpVyjAVcy5c8Dbb0vhw9tv5xe8kRFbQKGU\niKWz4blptmMW3fxmakQikpJmmrGPj0s0vWOHv/tgrwdSr9CCqFKGh4FbbpGUMCOAN9yQ+9F/+3bg\n+HHJ4Z2ZkfzbH/7QDr102gteHrAbYx1MT8sxzcJaY6O8trZKpExIHUIPuJZxD880WQx+Fr+GhmSy\nxIoVUhQxnwU+Lz/55EkbRQOS0UABJnUKPeBaZmREFsricclgaG62kaxXBZ1ZsLt4UcR30SKJnr28\n41xVeOa7DRvE6zWz30ZHgQ9/WKLfaFReN20q398HIdUI84CrlPZ2iTrPnJH82YYGmVRx773S5zcW\nk+/NwM2+Pms/XL4MvPcecNttsp17sc1r6OZzz2VGvRMTcs7mZuvffvGLc1PnCCHZoQBXAV7pZqZA\n4t13bWXawICIphn/PjQkvjAg+w4N2RzeZFL84uXLZRGvt9ce2wzdNNsZX9hkUZiJydPTIuAmc8FZ\ngFFoihwh9Qg94Crg4YdlioWJSJ1j3Xt7gSNHRDCNHxyJiAhPT8tCW3OzRKyjoxK1zszIdrGYfD85\nKccyx37xRYmqDcbL7e2VDmqmac/MjFgZPT1zBdbPbDpCahh6wLWCaZqTSsmriUgBiS7NZAqlRHwj\nERFHpWwXM+d7QEqLW1rETkgmM0fMb9rk7eUa2yMWk+NEIrLw5xXdstk6IfmhAFcJWssf04nMLI71\n94v/GolI2ldrq0SlixaJxwvYKLS7WyLglhaJiru6bA9gsx0gBRYPPiijirZtk/eAnOu220SUr7tO\nxhCdP+9dAMLiCkLyQwuiCjDz2cbGxDpobZVHevNY79U3whRJnD4tkW1TE7ByZebctulpSR0z0XEi\nAXzkI9beyIaf0uGgelkQUiVwJlwtYFLAzNy1REJeT56URbThYe/uY84uZs6c3z17pIXkxYuy+NbV\nJRMxkklgwQIb7ebCT6e2Qru/FQIX+EitwAi4wnEuZp06ZVPMTMWZu2FOIcy3CY5XQceiRSLGExPB\niyIX+EgVwEW4WsC5mNXZKT83NEiWQ1eXiKC7aMLvOKP5+rQmul22TOyMyUnJ0jh4sDwz77jAR2oF\nWhAVjvNxPxYDbrzRerhaSwQ6OZk57l3ruSPgvSJEY1PkK5zI9sjvFELTSQ0IXhTL0ayekHJAC6LC\ncS9m7dxpPdz2dumIduWK+MHJpHi+a9aIKM/M2OY8Dzwwf1sg2yO/2x4BpJuZ1vL6la+U7u/BCRf4\nSBXAZjy1Qq5FJ69G6q2tEiW//bYIcTwOrFo1f680m1fsFMJEQr4bH6coEgJmQVQnXmLrnLnmthT6\n+6UxjqmS6+yUHN3ly0WA43HrHWezBcw5z52T3hJdXSKeRkSzPfIHmelASD3ACLjM5Euh8nrcv3gx\nd7aCex/z+O/1uRF0L4E3EbPJGc6VZ8zolpCc0IKoRPKlUHk97re3ewusIZtA5irQ8BL4EydsJL1m\nDWezEVIEFOBKxC2w8bjNKPCaEGyi1lJFoFu22AnGkYiIbTQqC3Zai5dsPOMgF9IIqXEowJWIOwJ2\nlgZrLdVqiURxYpvL5ujstOljZpzRbbdJtkQ0antEOD1gQkjBUIArEbctcO6cbQcJ5H/s91OGm8vm\n2LpVIuBkUqLteBxYu9bfuQkhvmElXCViMgf275fH+6VLC6tG8zOKPlel2NKlssC2cqW8n5yUhj3T\n0+UraBgelgZDHR3y5+GH/Y+wJ6SWoACHTKEj2Z3iOjMjQm6E7KGHRMicJcbT03bixfbtUshx883S\nRjKRkNaUU1PAW29JNJ6rdLlU9PXZ7m5Xr0oZc5Cly4RUKrQgqgx39dkHH8jPZjLyxz8uHc2MzXHi\nhEy9iMdFaEdHJcPh6FFg8WJp6HPqlFgSt95anuY2vb3A4cO2fDkald7DtD9IDUELolLw2xzHD86I\n2XRFm52V11RKJlo4bY41a0R8AWncMzoq9oWzTDmVkmMB5WluY6yOyUn5MzEhi3+E1BsU4DLgx7fN\nhdMz3bBByn2ffVZmsZlRRGYwpxunHeEU2q4uiTzPnpXIOJWyqWhBe8H9/VIubUYomR7HhNQbFOAy\n4Ld9YrZIOZtn2t8vU4+VEjFtabHz2wzOiHnhQilRBuxI+c5OsR4iEfGFjQ9dyqjdzbJlkvp2113A\nnXcCt9wiv1QIqTfYC6IM+G2fmK3nw8iIHbwJyM8XL4qQHT6cu6Wke1S8c9umJol+GxtFeBcutIUX\nTq85V0vLoP9OCKllKMBlwG/f3WyRcnu7RLjJpI1259MQx72tO1/YKYJBNz33+3dCSC3DLIiAyFUw\nke27bE11hoeBHTvsOPpNm+zsNvP57Kz4u6tXS66vnyo2Z0Tc3CznNe0kvUqiWZZMiG9YCRcmuarR\nsn330kvAI4+Iz9vSAuzdC9x7b+5z7NtnswkAEeFbb5UR9IVYBu5rKkVJNCF1DPsBh0muR/hs3z35\nZGZfiKeeyi2ixhuembGfJZMS2XZ0FHe94+P5x9MTQoqDWRDzJF+WQK6Bl9m+8+O7Os974oQVa0Mk\nIiJc6KLWfAd0EkLmDwV4nuTL7c1VYpztOz8i6DxvW5tYFbGYbSsZj8txcw3Y9PrFUWhJNCGkeOgB\nz5Nsc9KKwc/kCa/zPvus/37B+RrCE0JKAj3gIClVHquf9pL5zltIKpqXzVHoNRBCSgMj4HlSqjlp\nhUakxZ7XK9UNYFRMSIkJPw1NKfVnAB4C8J7Wen2O7apOgEtFLisjiMjUS8AffbT0dgohdU5FWBBf\nAPDHAJ4N+DxVSy4rI9c4+vniZVewLJiQcAg0C0JrfQDApSDPUS3MJ/sgX1paqRrmVFIGRJBNgAip\nNAL3gJVSKwA8V2kWRLkXnuaTfZCtNLmYY1Y6tXhPpC6pCAvCN7t27br2c09PD3p6egI9XxCP97mY\nT3ObfA1rgm6YEwa1eE+k9hkYGMDAwMC197t37+7RWg9k3SFNRQpwOSj3/+jz8VndrSTdEXstere1\neE+k9nEHjbt27Rrws185KuEUfIbj5aTcpbfF+qyPPSaNdw4fltcdOyrLuy0VtXhPhGQj6DS0vwLQ\nA6AdwHsAntBaf8Fju1A84FLk8ZaLjg6ZiGEiw9ZW4MKFsK+KEJKF8POA/VItecBhVox1dEibSkNL\nS3ECzOo3QgKFAlxqnCv009PApUsydbgcAvbwwzILzowmammRuWrzPTezDQgJFI6lLzXOhbuhIake\nm++k40J55hngwQeB++4DFiwQ4S3m3Mw2ICR8KMAF4Fy4m5mR9o9AeQTMZETs3y9Rdzxe3LnZ/5eQ\n8KEAF4Bzhb6tzY54zyVgxVZ2ee1fCvFktgEh4UMPeJ74zaIo1mv12v/pp6srg4OQOqS6KuGqDb89\neIv1Wr32L6T/LyGkcql5CyLs5i7F2gX0agmpXWreggg73arYgo9qKxghhABgHrAQxOw2QgjJA/OA\nAT7CE0Iql6oVYL/eLtOtCCGVStVaEGF7u4QQkoPatiBYSksIqXaqVoDp7RZH2Ol5hJAqtiCYnlUc\ntHAICZTaroRjNVhx0MIhJHyq1oIg+cllM9DCISR8qtaCqEbKPYUil81AC4eQQGElXKVRbt+VVYCE\nhEZtp6EVSxhZAOX2XWkzEFLZ1K0A9/VJNHrxoox537AheCEutyCyCpCQyqZuLQjzeH76NDA1JeOF\nbrklWFuAvishdUNtp6EVS3s7MDoKJJMSjUajwdsCTJ0jhDipWwvCPJ43NcmAy85O+qSEkPJStxaE\ngbYAISQAmIZGCCEhwTQ0QgipZCjAhBASEhRgQggJCQowIYSEBAWYEEJCggJcBJwqQQgpBqahFQGn\nShBCssA0tKDhVAlCSDFQgIuA7R4JIcVAAXZRiK9bye0e6U8TUvnQA3ZRK75urdwHIVUKPeD5UCu+\nbq3cByG1DAXYRa34urVyH4TUMhRgF5Xs6xZCrdwHIbUMPWBCCCk99IAJIaSSoQATQkhIUIAJISQk\nKMCEEBISFGBCCAkJCjAhhIQEBZgQQkKCAkwIISFBASaEkJCgABNCSEhQgAkhJCQCF2Cl1MeVUseV\nUieUUr8V9PnCZmBgIOxLKCm1dD+1dC8A76eSUUr1+NkuUAFWSkUA/B8APwHgdgCfVUqtDfKcYVNL\n/xEBtXU/tXQvAO+nwunxs1HQEfBGACe11me01jMA/gbApwM+JyGEVAVBC/AyAIOO90PpzwghpO4J\ntB+wUurfAPgJrfWvpN//AoCNWuvPubZjM2BCSE2htc7bE7gh4GsYBtDleL88/VkGfi6UEEJqjaAt\niMMAblZKrVBKNQL4DIC/D/ichBBSFQQaAWutU0qp/wjg2xCx/zOt9RtBnpMQQqqFipgJRwgh9Uio\nlXBKqT9TSr2nlPpRmNdRCpRSy5VS31VKva6UOqqU+lz+vSoTpVRcKfWiUuqV9L08EfY1lQKlVEQp\ndUQpVfU2mFLqHaXUa+l/o0NhX08xKKUWKqX+Vin1Rvr/n01hX9N8UUqtSf+bHEm/Xs6lBaFGwEqp\nrQDGADyrtV4f2oWUAKXUYgCLtdavKqVaAbwM4NNa6+MhX9q8UEoltNbjSqkogIMAPqe1rvb/0fsA\n3ANggdb6U2FfTzEopU4DuEdrfSnsaykWpdQXAXxPa/0FpVQDgITW+oOQL6to0oVoQwA2aa0HvbYJ\nNQLWWh8AUPX/AQGA1vpdrfWr6Z/HALyBKs551lqPp3+MQ9YKqtqrUkotB/AJAJ8P+1pKhEIN9HJR\nSi0A8GGt9RcAQGudrAXxTfMAgFPZxBeogX/ASkQpdROAuwC8GO6VzJ/04/orAN4FsE9rfTjsayqS\nfgC/iSr/ReJAA9inlDqslPrlsC+mCFYCuKiU+kL6sf3/KqWaw76oEvFzAP461wYU4BKTth++CuDX\n05FwVaK1ntVab4Dkbm9SSnWHfU3zRSn1SQDvpZ9QVPpPtbNFa303JKr/tbSdV400ALgbwJ+k72cc\nwG+He0nFo5SKAfgUgL/NtR0FuISk/auvAvhLrfXXw76eUpB+HHwBwMfDvpYi2ALgU2nf9K8BfFQp\n9WzI11QUWuvz6df3AeyF9F2pRoYADGqtX0q//ypEkKudnwTwcvrfJyuVIMC1EpEAwJ8DOKa1/qOw\nL6QYlFLXK6UWpn9uBrANQFUuJgKA1vp3tNZdWutVkGKg72qtHw37uuaLUiqRftKCUqoFwIMA/iXc\nq5ofWuv3AAwqpdakP+oFcCzESyoVn0Ue+wEIvhQ5J0qpv4K0bWtXSp0F8IQx46sNpdQWAD8P4Gja\nO9UAfkdr/Y/hXtm8WALgL9KruBEAX9ZafyvkayKWGwHsTfdQaQDwJa31t0O+pmL4HIAvpR/bTwP4\ndyFfT1EopRKQBbhfybstCzEIISQcKsGCIISQuoQCTAghIUEBJoSQkKAAE0JISFCACSEkJCjAhBAS\nEhRgUjUopX5KKTVrkvbTk1bGlVIvK6WOKaV+qJT6RY/9vqaU+mfXZ7vSx1rl+Ow30p/dnX4fU0rt\nUUq9mT7+I0HfI6kvKMCkmvgMgH+CVBkZ3tJa36O17k5//xtOEU5X9N0NYEG6SZJBA/hReh/DzyCz\nouy/QHpI3Jo+/vdKeC+EUIBJdZAuud0C4D8gU4CvobV+B8DjAH7d8fFPQ+YQ/o3Hfl8H8On08VcB\nuAzgouP7fw/gvzuO/+Ni7oEQNxRgUi18GsA/aq3fgrQv3JBluyMA1jrefxbAX8FbgD+A9CG4HRIJ\n/435wvTCAPB7aYvjy0qpG0pwH4RcgwJMqoXPwgrklwH82yzbKaR7/iqlbgRwi9b6B1rrkwBmXG01\ndfqYn4EI/F7YxlANkFacB7TW9wD4IYD/VbrbISTkZjyE+EEptQjAxwDckW5AE4WI5594bH43ZBoJ\nAGwH0JZuQ6kAXAcR8t91bP9NAP8TwCGt9ZhSor9a6xGl1FWt9d70dn8LsSQIKRmMgEk18LOQuYEr\ntdartNYrALwNoBOOVqbpRbanAPzv9EefAfAT6X1WArgXLhtCaz0BYCeA3/c473NKqY+mf34AtdEm\nkVQQjIBJNfBzAP7A9dn/A/CfAaxSSr0MoBni6f6h1vovlVIrAHQ5B4lqrd9RSo0qpe6DYzSR1vor\njuM62wP+NoC/VEr1A3gfVd4mkVQebEdJCCEhQQuCEEJCggJMCCEhQQEmhJCQoAATQkhIUIAJISQk\nKMCEEBISFGBCCAmJ/w/VVy/4uVz8VQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df_data.rename(columns={ \"data1\": gene_expre1, \"data2\": gene_expre2 }, inplace=True)\n",
"sns.lmplot( x=gene_expre2, y=gene_expre1, data=df_data, fit_reg=False )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## BigQuery to Compute statistical association \n",
"\n",
"The Spearman correlation coefficient is defined as the Pearson correlation coefficient computed using ranked variables. BigQuery supports several statistical aggregate functions including a function to compute ['Pearson correlation'](https://cloud.google.com/bigquery/docs/reference/standard-sql/statistical_aggregate_functions#corr). Thus, the computation of the Spearman correlation is performed by the following steps:\n",
"\n",
"1. Combine the two tables containing the gene expression data\n",
"2. Compute ranked gene expression data\n",
"3. Use the BigQuery function CORR on the ranked data\n",
"\n",
"The following query performs steps 1 and 2. The two tables are combined by using an INNER JOIN command. There are several methods to rank numerical data, depending on what rank value is assigned to values that are equal, ['see this link'](https://www.geeksforgeeks.org/python-pandas-dataframe-rank/). Spearman correlation typically requires assigning the average of ranks to the similar values. Unfortunately, this method is not available as a simple BigQuery command. Nevertheles, the following query combines the RANK and COUNT BigQuery commands to generate the average ranked data. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"summ_table = \"\"\"summ_table AS (\n",
"SELECT\n",
" (RANK() OVER ( ORDER BY n1.data ASC)) + (COUNT(*) OVER ( PARTITION BY CAST(n1.data as STRING)) - 1)/2.0 as rnkdata1,\n",
" (RANK() OVER ( ORDER BY n2.data ASC)) + (COUNT(*) OVER ( PARTITION BY CAST(n2.data as STRING)) - 1)/2.0 as rnkdata2,\n",
" n1.ParticipantBarcode\n",
"FROM\n",
" table1 n1 \n",
"INNER JOIN table2 n2 \n",
"ON n1.ParticipantBarcode = n2.ParticipantBarcode\n",
")\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Once the ranked data is computed, we use the following query that uses the CORR command to compute the Spearman correlation coefficient."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" in runQuery ... \n",
" this query processed 7757877633 bytes \n",
" Approx. elpased time : 2627 miliseconds \n"
]
},
{
"data": {
"text/html": [
"\n",
"
\n",
" \n",
" \n",
" | \n",
" npatients | \n",
" correlation | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 540 | \n",
" 0.918271 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" npatients correlation\n",
"0 540 0.918271"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"query_corr = \"\"\"\n",
"SELECT \n",
" COUNT(*) as npatients,\n",
" CORR( rnkdata1 , rnkdata2 ) as correlation\n",
"FROM summ_table\n",
"\"\"\"\n",
"sql = ( sql_data + ',\\n' + summ_table + query_corr )\n",
"\n",
"df_corr = regulome.runQuery ( bqclient, sql, [] , dryRun=False )\n",
"df_corr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To test our implementation we can use the 'spearmanr' function available in python:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SpearmanrResult(correlation=0.91827140086976378, pvalue=1.2724813661443188e-218)\n"
]
}
],
"source": [
"print( stats.spearmanr(df_data[ gene_expre1 ],df_data[gene_expre2]) )"
]
}
],
"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
}