Pre_processing_Kumasaka_et_al_tables_step1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "43d4471e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.pyplot._IoffContext at 0x14c67d4a30a0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import gc\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "\n",
    "import pandas as pd\n",
    "import rpy2.rinterface as rinterface\n",
    "import rpy2.robjects as robjects\n",
    "\n",
    "import tqdm\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "\n",
    "from itertools import compress\n",
    "from Bio import motifs\n",
    "from Bio.Seq import Seq #, IUPAC\n",
    "from collections import Counter\n",
    "from os import listdir\n",
    "from os.path import join\n",
    "from rpy2.robjects.packages import importr\n",
    "from rpy2.robjects import pandas2ri\n",
    "from scipy.stats import ks_2samp\n",
    "from statistics import mean, median\n",
    "\n",
    "import session_info\n",
    "\n",
    "%load_ext rpy2.ipython\n",
    "%matplotlib inline\n",
    "pandas2ri.activate()\n",
    "plt.ioff()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "da30de1d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<details>\n",
       "<summary>Click to view session information</summary>\n",
       "<pre>\n",
       "-----\n",
       "Bio                 1.79\n",
       "matplotlib          3.5.1\n",
       "numpy               1.24.2\n",
       "pandas              1.5.3\n",
       "rpy2                3.4.5\n",
       "scipy               1.9.1\n",
       "seaborn             0.11.1\n",
       "session_info        1.0.0\n",
       "tqdm                4.64.1\n",
       "-----\n",
       "</pre>\n",
       "<details>\n",
       "<summary>Click to view modules imported as dependencies</summary>\n",
       "<pre>\n",
       "PIL                 8.2.0\n",
       "backcall            0.2.0\n",
       "beta_ufunc          NA\n",
       "binom_ufunc         NA\n",
       "cffi                1.14.5\n",
       "cloudpickle         2.2.0\n",
       "colorama            0.4.3\n",
       "cycler              0.10.0\n",
       "cython_runtime      NA\n",
       "dateutil            2.8.1\n",
       "decorator           4.4.2\n",
       "future_fstrings     NA\n",
       "google              NA\n",
       "hypergeom_ufunc     NA\n",
       "ipykernel           5.5.5\n",
       "ipython_genutils    0.2.0\n",
       "ipywidgets          7.6.3\n",
       "jedi                0.18.0\n",
       "jinja2              2.11.3\n",
       "kiwisolver          1.3.1\n",
       "markupsafe          2.0.1\n",
       "matplotlib_inline   NA\n",
       "mpl_toolkits        NA\n",
       "nbinom_ufunc        NA\n",
       "ncf_ufunc           NA\n",
       "numexpr             2.7.3\n",
       "packaging           20.9\n",
       "parso               0.8.2\n",
       "pexpect             4.8.0\n",
       "pickleshare         0.7.5\n",
       "prompt_toolkit      3.0.18\n",
       "ptyprocess          0.7.0\n",
       "pyarrow             4.0.1\n",
       "pyexpat             NA\n",
       "pygments            2.9.0\n",
       "pyparsing           2.4.7\n",
       "pytz                2021.1\n",
       "ruamel              NA\n",
       "simplejson          3.17.2\n",
       "six                 1.15.0\n",
       "sphinxcontrib       NA\n",
       "statsmodels         0.12.2\n",
       "storemagic          NA\n",
       "tornado             6.1\n",
       "traitlets           5.0.5\n",
       "typing_extensions   NA\n",
       "tzlocal             NA\n",
       "wcwidth             0.2.5\n",
       "zmq                 22.1.0\n",
       "zope                NA\n",
       "</pre>\n",
       "</details> <!-- seems like this ends pre, so might as well be explicit -->\n",
       "<pre>\n",
       "-----\n",
       "IPython             7.24.1\n",
       "jupyter_client      6.1.12\n",
       "jupyter_core        4.7.1\n",
       "notebook            6.4.0\n",
       "-----\n",
       "Python 3.9.5 (default, Jun 15 2022, 08:55:22) [GCC 8.5.0 20210514 (Red Hat 8.5.0-10)]\n",
       "Linux-4.18.0-425.3.1.el8.x86_64-x86_64-with-glibc2.28\n",
       "-----\n",
       "Session information updated at 2023-04-05 11:45\n",
       "</pre>\n",
       "</details>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "session_info.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfdd9eff",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "This script combines both python and R components and uses the Kumasaka et al dataset as a starting point. It first defines SNPs that are likely to be under cooperative control and also generates control ('no-effect') SNPs.\n",
    "\n",
    "1) identify a set of SNPs that are \n",
    "\n",
    "        * caQTL enhancers that have a SNP inside the enhancer that alters enhancer-wide accessibility \n",
    "        * no-effect enhancers where no there is a SNP that does not lead to large changes\n",
    "        \n",
    "2) once identified, we extract the two neighboring peaks or, in case the enhancer is within the 'dag' file, we will associates its 'dependents'.\n",
    "        \n",
    "        * identify all caQTL enhancers in the first column of the 'dag' file  (CM-leads) and extract their dependents 2nd column in 'dag' file\n",
    "        * add the two neighbors of all caQTL enhancers not in the 'dag' file and of all no-effect enhancers\n",
    "        \n",
    "\n",
    "\n",
    "Many functions in this script are sourced from the following script:\n",
    "**/home/.../custom_functions.R**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cfca7e1",
   "metadata": {},
   "source": [
    "## Switch to R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "e414a153",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "## genomic manipulation\n",
    "suppressPackageStartupMessages(library(Biostrings))\n",
    "suppressPackageStartupMessages(library(\"BSgenome.Hsapiens.UCSC.hg19\"))\n",
    "suppressPackageStartupMessages(library(stringi))\n",
    "suppressPackageStartupMessages(library(rtracklayer))\n",
    "suppressPackageStartupMessages(library(GenomicFeatures))\n",
    "suppressPackageStartupMessages(library(\"TxDb.Hsapiens.UCSC.hg19.knownGene\"))\n",
    "suppressPackageStartupMessages(library('AnnotationHub'))\n",
    "suppressPackageStartupMessages(library(org.Hs.eg.db))\n",
    "suppressPackageStartupMessages(library(GO.db))\n",
    "suppressPackageStartupMessages(library(ChIPseeker))\n",
    "suppressPackageStartupMessages(library(DESeq2))\n",
    "suppressPackageStartupMessages(library(Rsubread))\n",
    "suppressPackageStartupMessages(library(\"tximport\"))\n",
    "suppressPackageStartupMessages(library(\"readr\"))\n",
    "suppressPackageStartupMessages(library(biomaRt))\n",
    "suppressPackageStartupMessages(library(vcfR))\n",
    "## visualisation\n",
    "suppressPackageStartupMessages(library(ggplot2))\n",
    "suppressPackageStartupMessages(library(patchwork))\n",
    "suppressPackageStartupMessages(library(gplots))\n",
    "suppressPackageStartupMessages(library(RColorBrewer))\n",
    "suppressPackageStartupMessages(library(plotly))\n",
    "suppressPackageStartupMessages(library(reshape2))\n",
    "suppressPackageStartupMessages(library(\"wesanderson\",lib.loc=\"/home/kribelba/R_libs/\"))\n",
    "suppressPackageStartupMessages(library(ggpubr))\n",
    "suppressPackageStartupMessages(library(tidyverse))\n",
    "suppressPackageStartupMessages(library(dplyr))\n",
    "\n",
    "# custom functions\n",
    "makeTr = function(someColor, alpha=100) scales::alpha(someColor, alpha/100)\n",
    "\n",
    "source(\"/home/.../custom_functions.R\")\n",
    "\n",
    "## color palettes\n",
    "enh_cols = c(wes_palette(\"Moonrise1\")[3],wes_palette(\"Rushmore1\")[4],wes_palette(\"Rushmore1\")[3])\n",
    "class_cols=c(makeTr('indianred2',80),makeTr(wes_palette('Moonrise1')[2],80),makeTr(wes_palette(\"Cavalcanti1\")[3],80),makeTr(wes_palette('GrandBudapest2')[4],80))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "49abb584",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R version 4.1.0 (2021-05-18)\n",
      "Platform: x86_64-pc-linux-gnu (64-bit)\n",
      "Running under: Red Hat Enterprise Linux 8.7 (Ootpa)\n",
      "\n",
      "Matrix products: default\n",
      "BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so\n",
      "\n",
      "locale:\n",
      " [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              \n",
      " [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    \n",
      " [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   \n",
      " [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 \n",
      " [9] LC_ADDRESS=C               LC_TELEPHONE=C            \n",
      "[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       \n",
      "\n",
      "attached base packages:\n",
      "[1] stats4    tools     stats     graphics  grDevices utils     datasets \n",
      "[8] methods   base     \n",
      "\n",
      "other attached packages:\n",
      " [1] lubridate_1.9.2                        \n",
      " [2] forcats_1.0.0                          \n",
      " [3] stringr_1.5.0                          \n",
      " [4] dplyr_1.1.1                            \n",
      " [5] purrr_1.0.1                            \n",
      " [6] tidyr_1.3.0                            \n",
      " [7] tibble_3.2.1                           \n",
      " [8] tidyverse_2.0.0                        \n",
      " [9] biomaRt_2.48.3                         \n",
      "[10] readr_2.1.4                            \n",
      "[11] tximport_1.20.0                        \n",
      "[12] Rsubread_2.6.4                         \n",
      "[13] DESeq2_1.34.0                          \n",
      "[14] SummarizedExperiment_1.24.0            \n",
      "[15] MatrixGenerics_1.6.0                   \n",
      "[16] matrixStats_0.63.0                     \n",
      "[17] ggpubr_0.6.0                           \n",
      "[18] wesanderson_0.3.6                      \n",
      "[19] ChIPseeker_1.28.3                      \n",
      "[20] GO.db_3.14.0                           \n",
      "[21] org.Hs.eg.db_3.13.0                    \n",
      "[22] AnnotationHub_3.0.2                    \n",
      "[23] BiocFileCache_2.2.1                    \n",
      "[24] dbplyr_2.3.2                           \n",
      "[25] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2\n",
      "[26] GenomicFeatures_1.44.2                 \n",
      "[27] AnnotationDbi_1.56.2                   \n",
      "[28] Biobase_2.54.0                         \n",
      "[29] reshape2_1.4.4                         \n",
      "[30] plotly_4.10.1                          \n",
      "[31] RColorBrewer_1.1-3                     \n",
      "[32] gplots_3.1.3                           \n",
      "[33] patchwork_1.1.2                        \n",
      "[34] ggplot2_3.4.1                          \n",
      "[35] stringi_1.7.12                         \n",
      "[36] BSgenome.Hsapiens.UCSC.hg19_1.4.3      \n",
      "[37] BSgenome_1.60.0                        \n",
      "[38] rtracklayer_1.52.1                     \n",
      "[39] GenomicRanges_1.46.1                   \n",
      "[40] Biostrings_2.62.0                      \n",
      "[41] GenomeInfoDb_1.30.1                    \n",
      "[42] XVector_0.34.0                         \n",
      "[43] IRanges_2.28.0                         \n",
      "[44] S4Vectors_0.32.4                       \n",
      "[45] BiocGenerics_0.40.0                    \n",
      "\n",
      "loaded via a namespace (and not attached):\n",
      "  [1] utf8_1.2.3                    tidyselect_1.2.0             \n",
      "  [3] RSQLite_2.3.0                 htmlwidgets_1.6.2            \n",
      "  [5] grid_4.1.0                    BiocParallel_1.28.3          \n",
      "  [7] scatterpie_0.1.8              munsell_0.5.0                \n",
      "  [9] withr_2.5.0                   colorspace_2.1-0             \n",
      " [11] GOSemSim_2.18.1               filelock_1.0.2               \n",
      " [13] ggsignif_0.6.4                DOSE_3.18.3                  \n",
      " [15] GenomeInfoDbData_1.2.7        polyclip_1.10-4              \n",
      " [17] bit64_4.0.5                   farver_2.1.1                 \n",
      " [19] vctrs_0.6.1                   treeio_1.16.2                \n",
      " [21] generics_0.1.3                timechange_0.2.0             \n",
      " [23] R6_2.5.1                      graphlayouts_0.8.4           \n",
      " [25] locfit_1.5-9.7                bitops_1.0-7                 \n",
      " [27] cachem_1.0.7                  fgsea_1.18.0                 \n",
      " [29] gridGraphics_0.5-1            DelayedArray_0.20.0          \n",
      " [31] promises_1.2.0.1              BiocIO_1.2.0                 \n",
      " [33] scales_1.2.1                  ggraph_2.1.0                 \n",
      " [35] enrichplot_1.12.3             gtable_0.3.3                 \n",
      " [37] tidygraph_1.2.3               rlang_1.1.0                  \n",
      " [39] genefilter_1.76.0             splines_4.1.0                \n",
      " [41] rstatix_0.7.2                 lazyeval_0.2.2               \n",
      " [43] broom_1.0.4                   BiocManager_1.30.20          \n",
      " [45] yaml_2.3.7                    abind_1.4-5                  \n",
      " [47] backports_1.4.1               httpuv_1.6.9                 \n",
      " [49] qvalue_2.24.0                 ggplotify_0.1.0              \n",
      " [51] ellipsis_0.3.2                Rcpp_1.0.10                  \n",
      " [53] plyr_1.8.8                    progress_1.2.2               \n",
      " [55] zlibbioc_1.40.0               RCurl_1.98-1.12              \n",
      " [57] prettyunits_1.1.1             viridis_0.6.2                \n",
      " [59] cowplot_1.1.1                 ggrepel_0.9.3                \n",
      " [61] magrittr_2.0.3                data.table_1.14.8            \n",
      " [63] DO.db_2.9                     hms_1.1.3                    \n",
      " [65] mime_0.12                     xtable_1.8-4                 \n",
      " [67] XML_3.99-0.14                 gridExtra_2.3                \n",
      " [69] compiler_4.1.0                KernSmooth_2.23-20           \n",
      " [71] crayon_1.5.2                  shadowtext_0.1.2             \n",
      " [73] htmltools_0.5.5               ggfun_0.0.9                  \n",
      " [75] later_1.3.0                   tzdb_0.3.0                   \n",
      " [77] geneplotter_1.72.0            aplot_0.1.10                 \n",
      " [79] DBI_1.1.3                     tweenr_2.0.2                 \n",
      " [81] MASS_7.3-58.3                 rappdirs_0.3.3               \n",
      " [83] boot_1.3-28.1                 Matrix_1.5-3                 \n",
      " [85] car_3.1-2                     cli_3.6.1                    \n",
      " [87] parallel_4.1.0                igraph_1.4.1                 \n",
      " [89] pkgconfig_2.0.3               GenomicAlignments_1.28.0     \n",
      " [91] xml2_1.3.3                    ggtree_3.0.4                 \n",
      " [93] annotate_1.72.0               yulab.utils_0.0.6            \n",
      " [95] digest_0.6.31                 fastmatch_1.1-3              \n",
      " [97] tidytree_0.4.2                restfulr_0.0.15              \n",
      " [99] curl_5.0.0                    shiny_1.7.4                  \n",
      "[101] Rsamtools_2.10.0              gtools_3.9.4                 \n",
      "[103] rjson_0.2.21                  lifecycle_1.0.3              \n",
      "[105] nlme_3.1-162                  jsonlite_1.8.4               \n",
      "[107] carData_3.0-5                 viridisLite_0.4.1            \n",
      "[109] fansi_1.0.4                   pillar_1.9.0                 \n",
      "[111] lattice_0.20-45               KEGGREST_1.34.0              \n",
      "[113] fastmap_1.1.1                 httr_1.4.5                   \n",
      "[115] plotrix_3.8-2                 survival_3.5-5               \n",
      "[117] interactiveDisplayBase_1.30.0 glue_1.6.2                   \n",
      "[119] png_0.1-8                     BiocVersion_3.13.1           \n",
      "[121] bit_4.0.5                     ggforce_0.4.1                \n",
      "[123] blob_1.2.4                    caTools_1.18.2               \n",
      "[125] memoise_2.0.1                 ape_5.7-1                    \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0433d898",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "## load in Kumasaka et al tables\n",
    "## change path to local directory\n",
    "dag = read.table(\"/home/.../dag.tsv\")\n",
    "peaks = read.table(\"/home/.../peaks.tsv\",header=T)\n",
    "vars = read.table(\"/home/.../lead_caQTL_variants.tsv\",header=T)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "88b81139",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAMAAABKCk6nAAAC9FBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkbGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiosLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////d8OfNAAAYhklEQVR4nO2de3wURZ7AW8mtUbhlVVw5VtzVBV1hBdn1zoVbPPBOTle9011Yz9NDHYTI0wAKylOU8DAQCBBAEBCMvCESokGzSoxRITwF5CWPsATy4pHEZJJM/XMzVT2TybzSXV01VV39+34+dIbp9O9X+qVf1dX10xCgNJroBgB8AcGKA4IVBwQrDghWHBCsOCBYcUCw4oBgxQHBigOCFQcEKw4IVhwQrDggWHFAsOKAYMUBwYoDghUHBCsOCFYcEKw4IFhxnCz4oHY3Qt90SliEUOY1R/3fvtW2QmSjWON4wX/V/piPGu76z8C3ZYmvCWwTcxwv+GFtE0LLtI+bvn6hdYm4NjHHKYJ/eLzdz5+7iFDli7e37VuEUPlfbuq2xCu4m6Zpk9Adv/agRdotVy+01jaiIk2lXdghgi+3/4f/6qV1cXv+Rfttr2vaXfD8Tvtlj2u9gnPu1V45eEQbgJDnX7U3R2pPINTwk+6im8sQhwierQ1F6A/a1lytUwN6OuHdj7SuP6JnySF6M8rU3vD+zuGftE38abH3w12tPKLbyw6HCP4/DTN9jjYQ/32yNhWhHX7B6dq7vi8naVq672c/TaHraIcIflFzZXrZM0tLxn9P9gku8At+T5vp+3KcpqX6ft6v1QpsKmMcIjhVS0IoI3lvtnZvPXq2zaolWrdaNNwvOE8b5P2dvQk3t2l90vvhpvaim8sQhwi+fMu1//OYdmdVYzft3t7X/FNZRWut0wPX+AXXXt/be23VQ0t9XXsIoVLfJZcyOEQwOtrvxp8/ew6hi//boW2/7xDa2adtlwV+wWhwq9Nohtax9vLN2nI0R8sT3VqGOEVwbM5dPybwufHOvgJbwhwQjEm+sdr/catWKLIlrAHBmPKlxf6PO9eJbAhzQLDigGDFAcGKA4IVBwQrDghWHBCsOCBYcUCw4oBgxQHBigOCFQcEKw4IVhwQrDggWHFAsOKAYMUBwYoDghXHguDStYAErHfzEvzh0xmAeHoe5yY4nX5bgBkvgGC1AcGKw0TwlUivvINgKbAsuGZa5+u0hE6T60JXgGApsCx44BP5FfUVBf1fCF0BgqXAsuB2NfhHQ8fQFSBYCiwL7r4J//jsvtAVIFgKLAsuvK1r/xcHdPvF7tAVIJgr7gpjkzlZv4puyF2asvjThrDvQXBEPLXnKg/vrvQays9at2zN6ei/WTR31oYfw74tXeByfYw8G1wu1wEj+VjdB1fmhH4DggNU1Zcc+3bzyq3n9k5PchG+qk4hHw5H22i3b+380N20Yabv6y++xNteMpCbleAvEwMfv56OeXzm8eNXUc1xhy6/LThV7f1cnjlpzuwZo1wjZsxovtyof14ZLcJK/JuF+86UFu/bF1gW4a1S3sXbZmZnb8/JyY65nMS+J+tMLua/h5SW1qH6Uocs3SeWz0jd0ViTtXBRnqdm1fjxQ11p50vPT3cljfd9Dlsm+z9nZKSlpmWEL6fgtcsyV696PzMzsFyEt5qUhrfdtnv3gUMHd8dcJnPrqhzRn35bm1B/6NBXi4a5xnyByueSA+7iFb5lzlryty3oS1c0Vuk/F0cLnolXhz7sc+Nvc77DP8LP0OHw64tWWXDpB6tShriGBWztSm0m7039Z3LGpGh+F1/RP5RGy1EzN+J11A/ek/CaRlQ4y5X+dyNNtSz43BGd0BVqCr6QnTVz6JBQXZMjW1xQ+rdIX09JSfukHlXnfpj7xa4Y0156io9F3EWv1JtpsGXBb7dKaI8JXaGc4Kqs6WPC1BJGBT7ha+MM/fh7HHnSyafRb015Z/720tp9+d/XxLfV1g/R67pE/l4xwduSwrQ2sVY/BbtyTnsNz6r8cZHLNXiCr5CH58C6rNLSRnHtti646pnI36sj2HNs9WuDYuh1TWioXOg1OnT6Ke+F16li372r28gFUDyAi6zYNJ78YExUuSNdwybOmbkQXwk1SDoNPAiOTuORzFei7rVz1q6/IrqBRgDB0SiOcoszNGXe6gpZDsAtA4Ijc2hsJLmvbYjxcEBOQHAEyua8FMHuS1tFt4sGEBzOsXC5SVPWXxTdLDpAcCiVn4ZcNQ+ZtC326x9SA4Kbk928P2PI1M32uZ6KCAgOYk9yyDXVBdEtsg4IDnBxVDO7L6+4LLpFLADBfuY133sPiW4PI0Cwzqzmt0RfiW4PK0Cwj68mDA62OylLdIPYAYJR3bZmT3mHKLPzYhwvuDGt2bE5Jc7P47njbMENm4Y30ztWocrAOo4W/G1In1Wa6AZxwMGCPSnN9Y7+XnSLeOBcwRebPzF6t0p0g/jgWMHzm/VaZSlUtb05DhV8IvjsO+qE6OZwxJmCdwT1ahQKHNMaB5wo+FLQvdEESQdDMsOBggubDs+DwuYlUA7nCS4KOvnaeKSGURwnuKJp9y0U3ZZ44DDBV8YH/I6w+VgcgzhLcG7T4XmL6LbECUcJzgzoHazEcBwjOEhwY9CYqzLRjYkbzhG8qEnvoCLRjYkfjhEc9K7RaAfcHQVwiODyoQG9Q8+LbkxccYbg1U277wrRbYkzThD8Y9Cgup2iGxNv1BfcOK1J78gzolsTd5QXfCzoye8O0Y0RgOqC9wfdHKn/6CgCagsuD3ryOyysbIgjUFrwuia9ro2iGyMIlQUHXV2Nuiq6MaJQWPDEprPvMcFNEYi6gpuOzyOMTH2vKsoKrgv4XSK0HaJRVvBMXe9wJz1ZiICqgv1z6aeKbIQMKCr4uO73LYFtkAM1BZf4z7/V4togCUoK9u+/rpXCmiANjARHGOMkTvDJQO+kqfIVamJZ8IXnfz+2tHtC57ASbeIE+8vdvOTM3ufmWBb86FPr/3TLwsbUh0JXCBO80d9/pd6EGxRYFvzTS+hMYh2qujF0hSjBbv8BukRMfsmwLPjXWWiVdhgV3Rm6QpTg9frxGS6gMZYFb7q+/a3p9wy/M6yYrCDBHn3qjbNCssuH9avossIqtHNiWPlgUYL1KoJThSSXEA4FoqtPYAb+hb5V9JzWT8COefeoJTgUiM76M6ZTX/pWUePWJxWdIiC3nCjWkzVa34FtUbMqLjARfCXSTCYCBHte1/2uintqabEsuGZa5+u0hE6Tw3qN4i/Y4y+58LLqU+eYwLLggU/kV9RXFPR/IXRF3AW7A4UGDZXGdgiWBbcjEyw3dAxdEW/BDSP8fgvim1huLAvuvgn/+Oy+0BVxFuwJ+J0V17yyY1lw4W1d+784oNsvwt4Lia9gT+D4/E4808qP9avohtylKYs/bQj7Pr6C/Y+QXO/FM6sNUOQ++ILfb0Yck9oCNQTv8vvdEL+cNkEJwaf8fp36hlkMVBDsedkFPdDRUEFwtn8MljNmnzSHCoL9O3Ds/xSHooDg7foYO+igjIT9BRfo+29xfNLZDdsL3qL7HReXbPbD7oKP+e+QTsUjmw2xu2D/HJRwhxQFmwv+XPc7GN5CioLNBevzYA2GMVjRsLfgA7D/toS9Bes78FH+mWyLrQXvJn6TuSeyMbYWrL8IHPZqMtCEnQV/QPyO4p3H1thY8GUXnIFbxsaCX4cd2AD2FfyZvgMf4JvG7thWcI3uN4VrFvtjW8H79FEc0McRG9sKnkcE53FNogB2FewhfsfwzKEEdhWsv8ng4KncDWJTwfVJ2G8SvAjcEjYVvJPswPM4plAEmwomO7DLycUYDGJPwXrBjUn8MiiDPQXrQ91hJHTL2FLwEeJ3NLcECmFLwdOJ4HJuCRTCjoIbyXR2MB2lEewo+GOyA3/PK75S2FEwmW9lKHRyGMGGgs+RHXgRp/CKYUPBqXCJZQIbCh4P90gmsJ9gD+mmzOcTXTnsJ3gt9juMT3D1sJ9g0k05gU9w9bCdYP1B4UIuwRXEboLr9Te+oaqZQewmWJ+yYTKP2EpiN8FpRPAJHrGVxGaC64nfsRxCKwobwRURplDgIvg9IriIQ2hFsSz4aO/Tpx5olfBgWK1AHoI9et0rKDxpGMuCH5hQ/9j42to3HgldwUPwCuJ3AfvIymJZ8I116JdXEWpsF7qCg2DPECIYZpU1jmXB/dagpzch9Em30BUcBOvTFs5gHlhhLAsu7trjkVb/8VCHr0NXcBA8kgguYx5YYaxfRXvyFqcsygl/i5O94BLYgc3DoX6wDnvB+nD306zjKg2H+sHr+2I6/ht9qyJDJuV4nXVYtbFRT5ab7MA7GYdVHBvVD84lgmsZh1UcG9UPJvPavco4qurYp35wEdmBoTahOexTP3gCEVzKNqry2KZ+sD4vFhyhTWKb+sH58KSfCtvUD56G/Y5gGtMJ2OU+WJ8XazXLmI7ALoL1a+gLLGM6ArsInghHaDpsIria7MBQwN00NhH8BREMBWRNYxPBU+EITYk9BDeSHXgNu4iOwR6C98BYHVrsIXgB9gvlNyiwheB6Mlx2ObOADsKwYLfZyAwF74draGoMCa5b+GSntu0fGnfeTGSGgt8ngpnFcxJGBK/oPbOgDFXve//hVy8bj8xQMJmbEmb/psGI4PWBR0W7TdRIYCh4BBYMJexoMHYOzsLLlaYisxN8DEbb0WNM8H1PFqN9ve43FZmdYNKN9RqrcM7CmODG9Nv7d1hubvZPZoL1R8EzGYVzGAZvk9Z07NPV5IhzZoL1bqzvGIVzGMYE/7HXAfT5Pc+ZisxMMHljBaa2o8OY4GW+o3OduSnWWQm+SnbgpWyiOQ4jgj8N/PXYKeORWQn+Bt4ptIIRwXMeWnbEjTxns58abOJ5DivBy7DfkWyCOQ9Dh+ir0x68tePNvxtsqjOYkWAPmbtwLZNgDsTww4ayRpORGQneTY7QJs4NQDDSPy5cCA8aLGFMcG6vLj5MRWYkeAz2O4VJLCdiTHDHqd8d8WIqMhvBlWQHPsgiliMxJvgWihpFbATvJYLNXgAAfowJ/utC84bZCN4C/dDWMCa4p3bTXXfffbepyGwEj8aCP235F4HIGBN8kGAqMhPB5eQIvZ9BKIdiRPDN3j/m39tkIng7EXyVQSiHYkRwgv7HHEwEf4j9wtxn9EgueB4WXMAgklORXDA5QofNwQUYxojgVnv37vX92WsqMgvBekc0zOBPjxHB7fyYisxC8CR4J8kqUj9saCSFRt+xHMjBSC1Y76f83nIgByO14PfhCG0ZqQWTEe9QaNQKUguG4XbWkVlwHRFs+s1kIAiZBe+HfkrryCx4KZyCrSOz4HFY8GarYZyNxIIvkVPwSYthHA4jwRGe91gWvJoIDp+JGjABI8Ftwr+yLDgF+33DYhSnY1nw8EQfWmJi6ArLgskOvN1iFKdjWXC1q+f+kpIbSkpCV1gVfIEIDq96CZiBwSF6S5ePeByi58PMHCxgcQ4+//AgDoLJgNkUa0EAJhdZnvnPhH9pUbBeJynPUhBA3vrBh7DfJIp3ZoBgONQP3peB6fMofau87MKCUy3FALj0ZH2/FvPIY3Qt0pmFBWdbigHIWz9Yn8T/WysxACRv/eC/E8HnrMQAkLz1g/+G/cLcOpaRtX4wedKwxEoIwIes9YPfglMwG2StH0xOwZVWQgA+JK0fXALj7Rgh6YiONdgv1HO3jqSCX4PxdoyQU3AtOUJ/RR8B0JFT8CkiGMpwWEdOwbnY70T6AIAfOQWTYpQf0gcA/MgpmAx5L6IPAPiRUnAdTM3BDCkFf0IE0+cGAkgpmLz4/SZ9biCAlILHYsHr6HMDAWQUXEWO0DBBJQtkFLwVhryzQ0bB5FmwuTprQBRkFAxvnTFEQsFn4RTMEAkFk/F2o+kzA0FIKHgbjLdjiISCU7Hgo/SZgSDkE6y/03CWPjMQhHyCz8CTBpbIJ/gw9juWPjEQjHyCD2DBG+kTA8HIJ3gJFryPPjEQjHyCk7Hgb+gTA8FIJ7iaXGOV0ScGgpFO8HEYzcEU6QTnYL8wexIrpBO8GQvOos8LNEM6we9APxZTZBNcT07B5fR5gWbIJvgoEfwjfV6gGbIJ/hI6Ktkim+D3sOD19GmB5sgmeDa8lcQW2QSPx4Jh/jNmSCa4GsZjMUYywXugo5IxkgkmY97fps8KhCCZYLIDb6LPCoQgl+CTRHDsNgFmkEtwGvY7nD4pEIpUgj3DsOAF9EmBUKQSXEqO0AfpkwKhSCV4N/Y7rJE+KRCKVIIzsOB36HMCYUgleDgWvIM+JxCGTIL1l1ZO0+cEwmAgmEzLfin0a/OCd0A/JXssC95/z7V3bEaoNuw3zQueh/3ONbsZEAvLgv8w053XPp+J4ClY8C6zmwGxsCy4rfemZstv3AwE6zNUlprcDIiJZcG//cS76J/EQPAJ7Hck1BtlimXB2a17laPK+3tYF3wK7oI5YP0q+uwHVQi514UNhDQtmMwxC/VG2cKhQHR9BWawWcFzseAvTW4FxIZDgegNfTEd+5hsCplj9oDJrYDYyNOTpb8YHFamFrCEPAWi90M/Fg/kKRCdDi8G80CeAtFvYMHbzG0EtIQ0BaL1Ce5OmtoIaBFpCkR/SwRDPxZjpCkQPQf7nWhqG6BlpCkQ/TYWnGFqG6BlpLkPJt0cxfT5gIjIIrgS7oL5IItgco01hT4dEBlZBG+CZ4V8kEUw6cfKo08HREYSwY1kSHQFfTogMpII1qfHukyfDoiMJILz4CKaE5IIJm8lQa0k9kgiOAUL3kOfDYiCJIJHY8EX6LMBUZBDcDn2C9NjcUAOwWRI9Az6ZEA05BBM5vFPp08GREMOwathzDsv5BA8HyZf4YUcgkk3B0y+wgEpBOsD7mrpkwHRkEIwKYb1On0uICpSCCYX0fPocwFRkUIwGVG5gT4XEBUpBJOe6L30uYCoyCDY8wo87eeGDIJLsF+olcQFGQT/gAXPpk8FREcGwUVY8Gr6VEB0ZBC8BAvOp08FREcGwa9iwYX0qYDoSCBYn+Gukj4VEB0JBF+GEZUckUDwTpicgyMSCJ4NF9EckUDwm1jwZ/SZgBhIIHgcFlxOnwmIgXjBFdBRyRPxgg/CNRZPxAs+jwW/S58IiIV4wWSi91P0iYBYiBe8Egs+RJ8IiIV4weQi+hv6REAshAuuJx2VZfSJgFgIFww90XwRLvgS9ptGnweIiXDBZHaOLfR5gJgIF7wUC95EnweIiXDBU2FMNFeECybXWFCwkBeiq66cgTdH+SK66koB9guzzHJDdNWVFVjwKiO/CtAguurKbBgyyxfRVVfImOgzRn4VoEFw1ZUq7HcUFNPhhuCqK4XQE80ZDvWDdQwJXo/9vmU0C2AaDvWDt/wZ0zlt7drlGcvWxlyuzMh41ZW8bdeuCnRlFyw5LGex78lykwrgS+dcvlxVXXU55vJoUtIgl+uC2+1BHjcsOSwHcevJ+tDQ1JOHUl1zTxhKAlDBryfLmGCAM/x6skCwFPDryQLBUsCvJwsESwG/niwQLAX8erJAsBTwG9EBgqUABCsOCFYcEKw4/ATn3NvXEP94E2sSbRDxeuYRb+gV8X/vXed4CTZKb0dGfIx5Jd2BJ2m2AsGcIoJgxSOCYMUjgmDFI4JgxSM6SHBfR0Z84irriC9QTVQUB8FXIKLAiHEQDIgEBCsOCFYcEKw4IFhxQLDigGDFAcGKw13w1/e1e76GTajse3/W80AgIpvA1cMR04jn+t3U8yjTiJs7t33sIn1E3oIbbttY8/g0JqHOtv64YcY9/oiMArsSEcuInq6r697ozTJi2fVbK596nj4ib8G5XRDK78wk1Ae9EHJfW6lHZBN4S8/EQBtZRMzrgVD9aZYRC9sjtP739BF5C17aH6GK65iEqir3/g+8wx+RSeDzXfYnIpYRFz/69K8eP80y4tVb3z/7p7H0EXkLTnnR+09aY9XzvqX9Zn9EFoE9D39UkohYRkxptfHKuAdYRkTpWuItZfQReQtePMD7Ty6BzSQ7lwbcnR+IyCLw/EHIJ5hhxEX/jlBtqwqGEXf+5lDdgu70beQteEd3hAo6MQnl/uch7qaILAI/06bNDVqbAoYRc3yCE64wjDjRe5nvSSijjsj9KrpDXsOAyUxCretR68UfkVFg3x7MMKK7w47GCQ+yjLj9zj3VqZ3pI3K/D97VvePAsNkfqHhV83HJH5FNYJ9glhELu7Xrd5ZpxPm/+lmfw/QRoSdLcUCw4oBgxQHBigOCFQcEKw4IVhwQrDggWHFAsOKAYMUBwYoDghUHBCsOCFYcEKw4IFhxQLDigGDFAcGKA4IVBwQrjpMFL+uNGu/7SHQrOONkwZ5ea9KfFN0I3jhZMDrQ8fZi0W3gjaMFoz5GalzbG0cL/rxrh32i28AbJwuuu2fn8l6iG8EbJwue+hzy3L9CdCs442TBjgAEKw4IVhwQrDggWHFAsOKAYMUBwYoDghUHBCsOCFYcEKw4IFhxQLDigGDFAcGKA4IVBwQrDghWnP8HOzW0G9gdo6cAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "\n",
    "## 2 ) Defining the different SETS: CM, LOCALS and SNPs with no change and no CM\n",
    "## Kumasaka et el. defintions\n",
    "## 1) CMs using first column in dag file\n",
    "allCMLEAD_tab = get_peak_var_tab(peaks = peaks,vars,unique(dag[,1]))\n",
    "GR_allCMLEAD = make_GR_from_df(allCMLEAD_tab,c(\"Chr\",\"Pos_Left\",\"Pos_Right\"))\n",
    "\n",
    "## Get 2nd order depenent peaks \n",
    "allCMDEP_tab = get_peak_var_tab(peaks = peaks,vars,unique(dag[,2]))\n",
    "GR_allCMDEP = make_GR_from_df(allCMDEP_tab,c(\"Chr\",\"Pos_Left\",\"Pos_Right\"))\n",
    "\n",
    "\n",
    "# subsetting CM-Leads to generate data that is very pure, SNP in peak, P_Lead * p_caQTL >0.5 & only Single Nucleotide poly\n",
    "SNP_f = sapply(allCMLEAD_tab$Ref,function(x){nchar(x)==1}) & sapply(allCMLEAD_tab$Alt,function(x){nchar(x)==1})\n",
    "\n",
    "CMLEADS_SNPinpeak_2 = allCMLEAD_tab[(allCMLEAD_tab$P_Lead>0.7 &  allCMLEAD_tab$P_caQTL>0.7) & allCMLEAD_tab$Inside_Peak==1 & SNP_f,]\n",
    "GR_CMLEADS_SNPinpeak_2 = make_GR_from_df(CMLEADS_SNPinpeak_2,c(\"Chr\",\"Pos_Left\",\"Pos_Right\"))\n",
    "\n",
    "plot.ecdf(CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center,col=\"grey40\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "3da45e9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "# Gaf full LOCAL Lead tab & GR,  the neighbor peak exclusion was added only for the Jan -14 update table \n",
    "\n",
    "pot_local_peak_ids = vars$Peak[vars$P_Lead>0.7 & !(vars$Peak %in% c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5])) \n",
    "                               & !(vars$Peak %in% (c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5])+1)) &  !(vars$Peak %in% (c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5])-1))\n",
    "                                 & vars$Inside_Peak==1 & sapply(vars$Ref,function(x){nchar(x)==1}) & sapply(vars$Alt,function(x){nchar(x)==1})  ]\n",
    "\n",
    "potLOCAL_tab = get_peak_var_tab(peaks = peaks,vars,pot_local_peak_ids)\n",
    "\n",
    "LOCAL_tabf_2 = potLOCAL_tab[potLOCAL_tab$P_Lead >0.7 &  potLOCAL_tab$P_caQTL > 0.7 & potLOCAL_tab$Inside_Peak==1,]\n",
    "GR_local_2 = make_GR_from_df(LOCAL_tabf_2,c(\"Chr\",\"Pos_Left\",\"Pos_Right\"))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ef6b7031",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "#  non-causal peaks  the neighbor peak exclusion was added only for the Jan -14 update table \n",
    "pot_noncausal_peak_ids = vars$Peak[vars$P_Lead<0.5 & !(vars$Peak %in% c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5]))  & !(vars$Peak %in% LOCAL_tabf_2$Peak) \n",
    "                                   & !(vars$Peak %in% c(LOCAL_tabf_2$Peak+1)) & !(vars$Peak %in% c(LOCAL_tabf_2$Peak-1)) \n",
    "                                   & !(vars$Peak %in% (c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5])+1)) &  !(vars$Peak %in% (c(dag[,1],dag[,2],dag[,3],dag[,4],dag[,5])-1))\n",
    "                                    & vars$Inside_Peak==1 & sapply(vars$Ref,function(x){nchar(x)==1}) & sapply(vars$Alt,function(x){nchar(x)==1})  ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "f834a1ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 56791    25\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "## select of potential no-effect peaks\n",
    "NONLEAD_tab = get_peak_var_tab(peaks,vars,pot_noncausal_peak_ids)\n",
    "NONLEAD_tabf0502 = NONLEAD_tab[NONLEAD_tab$P_Lead < 0.5 & NONLEAD_tab$P_caQTL<0.2 & NONLEAD_tab$Inside_Peak==1,]\n",
    "\n",
    "print(dim(NONLEAD_tabf0502))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "63d518e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 56791    25\n",
      "[1] 56369\n",
      "[1] 7625   25\n",
      "[1] 7625\n",
      "[1] 2051   25\n",
      "[1] 2051\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "print(dim(NONLEAD_tabf0502))\n",
    "print(length(unique(NONLEAD_tabf0502$Peak)))\n",
    "\n",
    "print(dim(LOCAL_tabf_2))\n",
    "print(length(unique(LOCAL_tabf_2$Peak)))\n",
    "print(dim(CMLEADS_SNPinpeak_2))\n",
    "print(length(unique(CMLEADS_SNPinpeak_2$Peak)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "76084350",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nOzdeWAU5f0/8M/M3leSzU0OCCEkCAgiCHIoKKilh2KRatXWUsSqVFvr91urtbbybWmtbW3tQaktFu3dn7WItqUetEhRBLWIhHCG3Ncm2Xtmdo7n98fGbAiQc3d2N3m//spOZud5MkneefLMc3CMMQIAgNTDJ7sCAABwbghoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaBhfDh8+zHHcBRdcEH25f//+yspKk8m0ZcuW6JE//vGPPM8fP378fFf49re/nZWV1d3drUd1YXxDQMO49sQTTxw/fnzx4sUzZ84kIlVVv/71r69cuXLq1Knne8sdd9whSdLjjz+uYzVhnEJAw7gWbQh/8YtfXLx4MRE988wzR48e/cIXvjDAW3Jycm655ZYnn3yyra1Np1rCeIWAhvRz+vTpVatW5eXlFRQUfOYzn+no6Ige93q969evLysry8rKuuqqq959993o8a6urhtvvDEnJ2fOnDlvvPFG73XmzJnzj3/8g4iuv/76Rx99lIi++c1vVlRUXHXVVUS0ZcsWjuPy8/ODwWB7e7vT6eQ47vnnnyeiDRs2hEKhH/3oRzp/4TDuMIC04vP5CgsLTSbT9ddfv2TJEiKaOXNmJBLRNO3SSy8lolmzZi1ZsoTjuNzc3La2Nk3T5s2bR0RlZWVz587leZ6Ipk2bxhj7xz/+MXv2bCK6//7733///ZqaGiK66aabogVpmnbZZZcR0Te/+c377ruPiFatWhX9lKIoZrN5zpw5SbsLMD4goCHNfP/73yeie+65J/py0aJFRPTCCy+8/PLLRDR16lRFURhjN998s9Fo/NWvfvXiiy8S0YUXXigIAmPs05/+dG9AM8Y+9KEPEdFf//pXxtgf/vAHIvra177WW9aRI0fMZnNmZqbVas3IyGhsbOz9VFVVlcFg0DRNty8cxiF0cUCaOXToEBH9+Mc/5jiO47i9e/cSUXV19eHDh4loyZIlBoOBiH7729/KsvzZz3727bffJqJPfOITVquViKIBfU6dnZ1EVFZW1ntk2rRpDz30kM/nE0Xx29/+dnFxce+nysvLVVX1er2J+SoBiIiMya4AwPAYjUYiuuOOO6688sreg1VVVa+++ioR5eTk9Ds/EAj0fWm32893ZYfDQURdXV19D0YikX4fREXTfICrAYweWtCQZmbMmEFEJpPpxhtvvPHGG30+34EDBziOmz59OhG9/PLLiqIQ0W233eZyuX7zm99UVVUR0XPPPSdJEhH98Y9/PN+VJ02aREQnTpzoPXLw4MHvfve7OTk5Tqfz4Ycfrq2t7f3UiRMnCgsLLRZLor5OAMJDQkg3Pp8vLy+P5/lbbrnl2muvJaIpU6YEg0FVVS+66CIimj179rJlyziOmzBhgsfj6erqijaNp06dunDhQo7j6Dx90KIo2my2ZcuWRT+lKMrcuXOJ6Iknnnj44YeJ6Kqrrop+KjpupPdxIkCCIKAh/Rw7dmzlypVutzs/P//Tn/50U1NT9Hh7e/unPvWp4uLizMzMlStXHj58OHr89ddfX758eWZm5syZM3/2s5+dL6AZY3fddZfBYKirq2OMffe73yWiiRMniqLo8/minSdPP/00Y+yHP/whEe3atUvfrxvGHY4xlpSWO0AKam5urqiouOeeex577LHznaNpWmVl5eTJk6PjRgASB33QADFFRUUbNmx46qmnwuHw+c556aWXTp48+c1vflPPisH4hFEcAGd48MEHp02b1t3dfb4RGm63+89//vOCBQt0rhiMQ+jiAABIUejiAABIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUWm2aazH49m1a1eyawEA0IPn+WuvvdZkMiXi4mkW0K+99tr27duXLl2a7IoAABARPfPMMxdddNGUKVMScfE0C2giWrx48R133JHsWgAAEBG99dZbibs4+qABAFIUAhoAIEXpHdCBQIAxpnOhAADpSI+AFgRh06ZNVVVVVqs1IyPDbDZXVlZu3LgxEonoUDoAQJrSI6A3bNiwf//+rVu3trS0yLLc3t6+bdu26urqu+++W4fSAQDSlB6jOHbs2FFfX2+z2aIv3W73woUL58+fX15erkPpAABpSo8WdGlp6c6dO/sd3L17d25urg6lAwCkKT1a0Js3b16zZs0jjzwyffp0l8sVDAZramo8Hs/27dt1KB0gdTDGOI6Ly6VEUbJaLQOc4PMGHE6b0Ri333FVVTvau3Jy3SZT+s2fSFN63OgFCxbU1tbu2rWrrq6us7PT7XavX79+6dKlBoNBh9JhPGhqbLVYzI31rW/ufcfb6TMa+LnzZ3k83Xv+9ZYkSBFZVhW1ZGLRDTd/+N0Dh5saWqdOm3z9J1a2tXb8+bcvdnf5Zs254PobV7a1dPx+219rTzaEAiGLzVJWXsIR11Df4nDYPvmZVRdfcmFvcaGI5DCfEY6Kovxg0y9e/tvrPq9/SuWkRx+7v/KCKUTk8wZ+8v2nt/95Z25+tqe9K3ryvf/72VvXrR7xF/vu/vc3/+iZ/x44TET3Pbj+k7et6nfCkfeP//CxX767/30iuvFT1973lfW8YbT/K7/4/CsbH3wi+vE5C4VE4JI16M3r9e7bt++aa64Z1rv+9Kc/eTwePF0cPyQxIkrSnl37Xv7b612dXuL40ydOi1KEGBGRYuUjNgPPNKtXFTONwSJTxGo0h5TskyKvxn6wGdHZrdYZs6sOHzwqZRjap9k5RitsZTW7D52vGp5KW/t0x8oLLnr4qtVf3vGb/9QeXXnBRc/e8nmrsWcFho0PPfHiX17p+5a/vf5sbl72l+58dM+/zjHT7Ikt31i89JIR3JCuTu+HFt/S98jjP3l46YqFvS8lMXLZRdf3+Tz7wgO337L24yMoq9fxmtpbVn2+75Gf/XrTvEtnj+aaY8btt9/+4IMPJmiqd9ICeu/evcuXLxcE4XwnvPXWW2evi/TGG2/cdNNNhYWFRDR58uRJkya1trbW1NREP4sjaX2EMXbiSNN//v2OFBYuXTar58fg9feO19SXlhUuvPriEK/6OeWdmpN1zW00K39hQWn0nP2dTcHjHktF9qLcniOnjzXK/zw9sazwsqvm9b3O2UeED5fOKp9ERE5m9OypPXmisbQ0f8nlc3rO2Xf45InG/KkFeUt6Hmg3iMF3PI15Vsei3NLJGTnZFpufZ8xsfvWF15bMmt7zriPHTze3Lb1s3gUTCrs6vb1HivKye8851tpWWF5iISq02qNH/LJ8ZdH0SESMqD2/FEbeYrXY+x05XlO/Z9cbfb8Kq9W+8vplvEGNHvG0+X754z/1/UprjzcvvXKxymTeoEZ/3RWZjLxJ01TepBExIpIlxnMmRqrB3BMISoSRZmCcYrLwwUDI1x3oew8dTrsrwxH0iRFRMZmNLnfP/xMBr6TKGmegzGwbY8Rx5O8WNIWMJt6RaYn27gR9kiTIFqvJmdXzLn+3IEua0cxlZtuJiDEK+kRZ0kxmgzPLwnE9R+YvmmO1mRsaGqI/LUVFE/Lz87u6uqNHiCg3N7e4uKjvkQkTCvudE/cj7777zrx5l4y1tTgWLVo0QDoT0YQJE+bOndvv4FtvvdXY2Dh9+nQiysjIICKn05mdnR39LI6k1xG/N/Sv42/s/NsuTVXkiMyIHas+3eLzWYoy8k818honmtiBSXLz1JzTFq45Uhd9V12u7M902FnkVLA7eiSoRPwlFrsSO9LoiEQWZXZaWTBwOnqktlT2Zp/jCLMrve8KucTuqTbJJhu6en4JT7nE7qm2kEsOfHBOVzhERCElcirYnRkxGHjJGxTCYUltDbVntkfPUVqCNq8iNHSHzK72uvbeI6oaO4cPS0W8g+c44wct/Syy+v1BIk3RVEYa00hVIp2yj0g1mnmNMSImCf7W1lafL3T6RFP0Xd7uQJab93Z5zTYDMSKO83n9RNT3nGBACAR8BiNnMHJExHGcKjNZk40mgyZrHMcRcYyYqkYMBl6VVSLiODIZzArTOI4njQt0h06faPJ2B3qvXDalJL8gx+nI4IhjHCmK1PNdLskyGs0RJcKYrKqageczXdlms1lVVUEK8zxHjOW4czmO5zg+LAZ7vvaMbKPRqKpaRBGJiBjLzc41Gs2apoWFAHEcEcvLLcjPzzWaDMFgkIg0TXO73TabLTubokeIKCcnu9+Rs8+J+5HOzp6eq0RIWgt6ZO67776WlpY//OEPya4IjFA4LGx6+Ecv/+M/nqlWb5HZ5lVyT4rtVXZmoJzj4dYLXWJWz5OJqTu7Tl+WKduH9KDC5lWErFhrI786lHtsoD//vTqm2Tum2Qd9S/Nsh3ey7dJJU79/3ad+9p9//vbtPfdetvJbH76x94nfpz5+79Hqk33fsu3//XDajIo7bvnywXeqz77g/V+7/eL5MxlTNU2LHon+IkZbiwYDzxgRMZ7nOY632WyKohqNPM8bwiHxkx+9p++lHv/p1y6/ckFvTWRZ2fCZh/779uHeE776f/det2Z4fYn9tLZ0XHvFZ/oe+d32n1RUTR7NNceMhHZx4GksJNyf/rHr8a3PHC9WNAPn6JAzmqSma3sa0VKmMVhoViw8EQUKzX3f1XhpxuDpzMgcUuweJe+YINv5YK7JoDJbp2zvUgat1dUfvry9vfO/Bw5n1YmM6PO333yk+Pjru/b1uXas5/qmnJmPPvK/DruNiLasWb9lzfp+V/u/73353nUPt7Z0RF9ef9M1vEk5evTImk9dIwjCsSO1fU/+yPVXLLp8ntHIWyxWm81usVh4fhgP8Z763eNbnvzNgTcPEtH/fu2upcsv7ftZk8n4tU1f/NFjv9z92j4iuud/1o4ynYmocELeU797/Fc/+/2be965ZOFFn1q3GumsDz1a0M3NzYFA4JyfqqqqGtal0IJOC5IiP7bzrz/d+WLYyojnRvYTVvROsLPCKmX0tCE4jZnCzCCrjg4lq04UXUZzWLH61b5v4XiOI46ImUym8oqJDqej/nQjEZ+R5bDZbVOrJn9q3eoTx043NbRWTps879LZiqLs3f12l6f7wjkXTJk6iTF28O3q5qY2p9POGCubUpqZlXHsyKns3KyKyrJoEYyxUCgUCgUEQVBVVdMUVSUijTFiTDt+pM7vDc27dFZ+Yb7FYrZYzCaThTHW3Nhmd9isNounvatk4oS4jLSTZQXD3VJB2regn3nmmYcffpjjuLNnprS0tOhQAdBHfXfnXT/8wV5fg2LlGBGzE3eO0RP9x1TYuhQh20hEOSfCYpY5lGvkNMqqFbLqxYwm0V9q1Xgus1E0RBhHZDAabU6rypmdzH79bR+aPqsyPz+3vHKSJEastoEGBfcqLi3s/dhoNF5+5YLelxzHXTRvxkXzZvQeUVW1akapzxc6caJGUdQP6s4ZjUaj0Wiz2S0Wq91utdkc0fOnT48Nxet72d5CSycVDaWSQ4F0Hg/0+B5/5StfqaioePTRRw8dOu8wJkhfQUlc88Nv/6ejVjNyRETWnvg9dyuRUfG7waaLnRwRr7Hsk2L+4ZBi4RlPJkEjChNHxMjIG8xO+9Jl89fefVNZeemgdRhiOg9KFAWfrysQEFRVYUzjOEbEG41mpzMzJ8dtMlnjUgrAEOn0R3jlypWYNzj2PPnYr7737j+7yqxERMYzArlfM9kQ0SxB1RzUCt4PGiIsq17se7IpQiaToaSi5MZbr73+ppU61Lyvzs5Or7c7urYizxPPG6xWW3Z2gcPh0rkmAP3oFNAOh+PZZ5/VpyxIKDmiXHvlZzo93YEic+PFLlZ2jkZlNJ0NMstslNwnBUtQPfscjqhgQv4v//B4foHeS7KoqtLe3u73+zVN4zjGcUaHw1ZaWmo2x6cZDhAv6MaCoXr3wKE7P/0QY1r7NLtncS47ZwczkTmoutrkvMMhXuv/dJDjuMKi3NvvvOVja67SocL9BAKBjo62SCTCGBmNxqysrJyc3DguVQEQd/jphMFVv3/s9hv/RzKyxkWuUG7P5GauXz+GxoreDWY2Sv0GbXAcd8WKS7/27f9xOJPTgStJUnNzgyRFOI63263FxSUWC7qSIT0goGEgiqJesWCNFJJql2UKWabe49EQ7k1nq18uf83LolO+iDgih81+xxdvSeKSOoIgtLa2RiJhxsjhsJeUTDWZTIO/DSCVIKDhvI4frb3lus9rRu7Yx3I0wxn9GT0vGNk7IyX7g0ZJI+I4oqysrB9v3Vg1PSFjQodCEMJtba2CIPI8ZzZby8rKLRZbsioDMEoIaDi3qxfe4u32tl/g8FSdI+A4jdynhcL3QkTEOCrIz/76d/933qWzdK9mTFtbq9fbzRizWCxlZWU2mz2JlQGICwQ09LfjuZf/76s/JKKuCtvZ6cwrWsHBkLtB4oj77Iabb7/7puSu6x0Oh1tbWyRJNBoNRUXFLldGEisDEF8IaIhRVW35/JvCoRARaSZqm+Hod4L7tDjhv0Gb1fqPd56z2ZP8qM3v97e2tjCmWiy2qqppPI/9H2CsQUBDj+NHT91yXc8yafWLM4J5ZuozToPTaMour0NkL+79rTs7K3nVJCLq6Gjv6vIQUXZ2bl5efnIrA5A4CGggIjp0sGbdjfcTkeA2nb48g3FnTNfmGF3wgudLX/3cTZ+6Nnl1JCJqbW32ev08TyUlkxyO/g18gDEGAQ10/TXrm+qaiSiYb65f5Dp7FY38Gmnn3t+7s5PYvcsaGupDoZDBYCwqmpCRkZm8mgDoBwE9rmmqtmT2KkVRiSiUb65fFIvgaOeGUdJ+cNmNn/3OR5NXR/J42j0ej8lkmjx5isWC2dgwjiCgxy/G2JJZq6LrGRNRw4IzGsgcUW595OQzvzEYk/bwzev1trW1cBxNnDjJbkeHBow7COhxStPYkguvVVSNiMQMY/2SDK1PDnNEP5n+sdu+c0OyqqcoSm3tKU2TCwqKsrLcyaoGQHIhoMej/W+8u2Htw9GPOyvsbTP7z+l4/qo7VyxfqHu9ejQ3N/j9AZcro7i4JFl1AEgFCOhx5xsPfP9v21+LfnxiuTvi6t+D8dDl1yYrnSVJOn261mjkKyoqjEbz4G8AGNMQ0OPLjR+++9SpOo4olGOqX5LJzlow9GtXffyB5dclo2rU2enxeDrc7pz8fAxtBiBCQI8r1191e1NDC0fUMdXRMePsFTa4+kd+km136l8xTVNPn66VZXnixElYQwOgFwJ6vLjpI3c2NbQQkeI0nJ3ODpFv+sEvjcmYLS2KYl1drc1mLy+v0L90gFTGJ7sCoIeNDz5x6mQDETGOji3vPyhiZc7Uth8+nZR0DoVC9fWncnJyJ06cpH/pACkOLeix79TJuheffzk6P/DoR3Kj8wR7F9lo//oWuy05yx41NTUEAsGiomLMDAQ4JwT02HfLdfdE9zo5dpVb++AbHk3nIw98L0npzGprT6mqVl4+xWzGaA2Ac0NAj3Ff+vy3VEXliDqr7IrjjE6MFz775VJ3nv5VUhTl9Olao9FYUZG0jVcA0gICeizbt+ftPa/sJaJAkbntAjv16dn47IJlV1bO0L9KghCur6/LyHBNmIBJKACDQECPWaFQ+J7bHyFimtHQML9nnY1oOl9SOuXJ69fqX6Xu7u729tb8/EK3G7O3AQaHgB6zViy4kYiIuNqlZzyCs3KGXRse0b8+Pp+3o6OltBTLHgEMFQJ6bPrEh+9UFY2ImIGkvpO5GbVt+qX+9WltbfX5ukpKkM4Aw4CAHoO+9bUnT59qICLG05GP5fb91L57Nxp4vQe/t7W1+f3eKVMqjUb8vAEMA35hxppD79Zs//NOImIc1ZyZzm6rfUax3vNB2tpa/H7f5MlTkM4Aw4XfmTFF07Tbb7k/+nHTPFe/tZCe/+z/6FyflpamUCgweXIF0hlgBPBrM6Z8aMmtTCMi8k20+ovP2B3qWytvnDdR13HHHo8nGAxVVFRx3FmL5gHAECCgx46vP/A9b5ePiEK5pqaLz1iUbtsnN6yePV/PynR0tHd3d06aVI50BhgxBPQYIYnS37fvIiKNp7rFZ4yrWzXjEp3TubOzs7u7s7y8wmg06VkuwBiD1ezGiLU39vQvn74si/q0WU0c/5tPfV7PmgSDgc5OT1nZFKQzwCghoMeCk8frTxw9RUQtcxyi+4z/ig59+Xt61iQcDjc3NxYXl2AJJIDRQ0CnPVVRb/7Y3UTkK7Z0TzpjJf637ttU4s7RrSaapjU1NRQUFDocmI0CEAcI6LR35SVrGLHgBHPTJa6+xxeWlE8vKNazJnV1tU6nKzMT62wAxAcCOr398se/FQRJdvL1CzL6HudVevnzX9ezJs3NzUT8hAlFehYKMLYhoNOYprGnfvp7Ijp1ef9Ga8c3f6FnTXw+bzgcmjQJ21YBxBMCOo194Y5vMGInl2ep5jPGGv/qkhstFsv53hV3iiK3tbWXlJTwuq/yATC24TcqXamKuu8/B+oWuc5YrI7otlDZjas/rGdNGhoa3O4sq7X/TuEAMEqYqJKu/u/BJ0Jucyj/jJZy/vHwT597VM9qtLW1chyfl5evZ6EA4wRa0OnqHy/+q36Ji4j1HrF2q3/+8tf0rIMgCIGAv7S0VM9CAcYPBHRaev6Pfw/mGRnPfbCJFZlEdcYBYe6CWbrVgTG1sbExNzffYDAMfjYADB8COi1995ub2yrtsdeMKnZ2//RXm/SsQ1NTi8Nhz8rK0rNQgHEFAZ1+3tj9tiqrYnbs+UH2ScFkMMy8qEq3OgSDAUmSMOoZIKEQ0OnngS98y1dqZnxsaF3+MeFzX7hVtwooitLa2jxhQiGWEgVIKAR0mnl6yx9EQWqZE5vVzSuaQWa3rf+EbnVobm5wOjOx/StAoiGg04kkSpufeNZXatN4rnf0RmZT5IabP6JbHTo7O1VVKyws1K1EgHErOQHd2dmZlHLT3Z2ffoCINc11sN7RG4wm/Df4Pw/fqU8FVFX1eruKi0v0KQ5gnNMjoNvb29etWzd//vyvfOUrHo/n4osvLiwsrKqqqqmp0aH0MSMSkQ+/d7xxfiZRbEX+zDrpqpWX69YX3Nzc5HJlms36zSMHGM/0COh169b5/f4HHnigurp6+vTpd9xxhyRJd91117333qtD6WPG59d+lYj5C2LblHCMSg6FvvWDB/SpQCQiaZqan49JgwA60WOq9+7du+vr6zMzM+fPn19ZWfnZz36W5/n169dv3LhRh9LHjIPvVEsuExlijeWSff7f/fXHulWgqanR7cZazwD60aMFnZ+f//rrrxPRv//9b1EUT506RUTHjh3Lzs7WofSx4Qff/gVjzFPVZ3KKQiWiZcrUMn0q0NXVZTQas7LwLQPQjx4t6Mcff/wTn/hEZmYmY+ynP/3p6tWrr7rqqhdffPFLX/qSDqWPDc/9/m9EFM6Nfb8y26Rf/78f6lYBv99XXKzr/iwAoEdAr1q1qqGh4cSJEzNnznQ4HLNmzXrllVd++tOfXnPNNTqUPgb4fAE5IjOeZGvsP57ik0pxiU5j3TyeDrPZaDJhH1gAXem03GhOTk5OTs/upUuWLFmyZInX6925c+cAGR0Oh1tbW/sd9Pl8jLFznj+Gffr6e4moY1psYgin0cYH79andE3T/H7fpEmT9SkOAHolbT3o6urqVatWCYJwvhNee+21Z555pt/B//73v+NtX6VwWGxpbiei7rLY4DZzSPvQx67UpwIeT4fT6cKSdQD6S1pAL1q0aIB0JqKPfvSjH/3oR/sdvO+++1paWhJZr5TzzYeeICLJZVTNfO/8lIWCTg/rVFUJBoPj7Y8iQIrQeyZhIBAYh30Uo/HqP/cSUcNCF32QzhyjZ76r07YpnZ2dmZkugwE77wAkgR4BLQjCpk2bqqqqrFZrRkaG2WyurKzcuHFjJBLRofS09qPvbmWaFnEaI/ZYD0NmgLJz9ViFWdO0UCjodufqUBYAnE2PgN6wYcP+/fu3bt3a0tIiy3J7e/u2bduqq6vvvlunx1zp63e/fp6ITl0RjeOe/zxeuushfUrv6up0OjOwVzdAsujxr+uOHTvq6+tttp5dn91u98KFC+fPn19eXq5D6enrtZ17mKZ1l1m1ntYzR0QmUZt9oR4L82uaisEbAMmlR+OotLR0586d/Q7u3r07Nxf/Ow/ka//7PcZRy0XOvgd/cdmn9Cnd4/E4HBi8AZBMerSgN2/evGbNmkceeWT69OkulysYDNbU1Hg8nu3bt+tQepo6+Ha1HJGb5kcX5u8ZvpHZGFnznat1KF3TlFAoOGlSmQ5lAcD56BHQCxYsqK2t3bVrV11dXWdnp9vtXr9+/dKlS9E6G8A3Hvi+ZuL8RdGxzxwRcYx9df6H9Sm9o8PjdLp4Ht8ggGTSafiUwWBYsWKFPmWNAZqqNTW2eqadsadU8bvBux/To3+DMSYIwsSJGPsMkGR4QJ+Kfv7ks0TUWWHtHTFuiGg3zFygT+nd3d02mx2DNwCSDhMQUtHvnn6+fbqDGWMbpRS9E3rkn3os/qdpmtfbVVqK5jNA8qGVlHK6unySLHdW2j44wIhRsWQxm00DvS1OfD6vw+E0mfQoCwAGhoBOOfes/Wr3JCuj3hnxXEaT+JVv6DGphzHm8/mzs3N0KAsABoWATjknjp32lVl7N4blNJr4TmjFNZfpUHQg4DebzWg+A6QIBHRq2fH8y4wx0WWgDxK68FBwzrwL9Snd5/Ni10GA1IGHhKnlD7/+q5RhYJqXLiIAACAASURBVH12hs2qk77zBz0W35AkieP43hn5AJB0aEGnltMnG08vzux9aZBZhtOemenSoeiurk6XS4+CAGCIENAppL3N01ZmVC2xb4qjPbLuzk/qULSqqrIcyczMHPxUANALAjqF3P+5jW1nzh4sPBi6/iY9pnd7vV6bzd77ZBIAUgECOlUwxt5vqGPGWES62mQbZ7TZLQO8K16Fh0IBt1unbbQAYIgQ0KnitX/u8ZaY+x4pfcM375KZOhQdCoWMRpPRiCfGAKkFAZ0qNv/wWd/EWGPZ2SYT0cObvqhD0T6fNytLjz20AGBYENCpoqmuRbbFlvfMPRF2ZTjz8hM+qU/TNE3T7HbH4KcCgL4Q0CnhWE1tIJtnvd8Nxuwd8s1rr9ehaL/fZ7Vi7DNAKkJAp4Stm3/fNdne+9IoMSK6bd0NOhTt9/sxug4gNSGgU8K+ve8GC829qz9n1Yt2l8NoTvhTO0EQDAaDyWQe/FQA0B0COvkYY10sohlig5Ddp8XLll2iQ9Hd3V1oPgOkLAR08u147uWW2bH+DY4xU1j73D23JrpcRVEURXU6nYOfCgDJgIBOvu3/b6eYE1vh09kSMRgNJRMnJLpcr9frdDowexAgZSGgk++Ap07r832YcChcUJinQ7mSJGD4M0AqQ0AnWXNjW/NkS28zlleYUVA/+enrEl2uIIQ5jud5w+CnAkCSIKCT7KH7Husd/syIMhskIrr6I0sTXa7X68vMRPMZIKUhoJPsSPUJxdHzXeCICg6HjCaDOyexIys0TVOUiMOB2YMAKQ0BnUxtLR2ymVPMHwS0yniFzV0wO9HlhsMhq9U++HkAkFQI6GT6wbef6p4UG79hjDAi2vSDBxJdbjAYzMjA5ikAqQ4BnUx7du3zlVh7X9q6FavV4spI7MBkTVNlOWKx6LDMNACMCgI6aTrau2RZUWyxb4H7tHjxJTMSXa7X67VasXkKQBrAGu1J86PHfhXKNWl9NvB2tEe++Os7El2uKIr5+QWJLgUARg8t6KTZvWtfR1VsnU9DRDMajWXlpQktVJZlxhg2TwFICwjo5AiFBDEsRFyxJ4QZDdLc+RcmutxAwOdwYPwGQHpAQCfHM0/9P+JIiT0gpIKa8D0PrEt0ueFw2OXC8nUA6QEBnRz/fnmv5DD2PqkzyMyk8ZVVkxNaaCQimUxmgwHTuwHSAwI6OZqa2oITYv0b5pBaOS2x6UxEgUDAZsPuVgBpAwGdBKqqSqLkL4rtY2L1ybfd8YlElyuKYacT81MA0gYCOgn2vv42EUUcsaEUGY3y5csXJLRQURR53sjz+I4DpA38uibB37e/ppk41RwbAe0OGxI99C0YDLhcaD4DpBMEdBK8tffdYE6sf8MQYROL8xNaImNMFAW7HcvXAaQTBLTeVFX1+4JGRes94mqWLlk4J6GFiqJoMpnRvwGQXvAbq7cdf3mZiPyFsRZ0ZqN0w80fTWihwWAQjwcB0g4CWm//fvVNIvJNii0mp7iMpZMSu0WsKAp2OwbYAaQZBLTeTh2rYxyppp47z4jKuYyEligIgsVi4Th8rwHSDH5p9dbV2e0vik3xNsp0202rElqi3++z27H+BkD6QUDrqrvLK0lyZ2UsoE2CsnjpvMSVyBiTZRkd0ADpCAGtq5/8YBsRMUPstmc0REonFSWuREEQTCbz4OcBQOpBQOvq1X/sYTzJtlgH9MSWxJYYCmF+CkC6SkJAd3d3BwIB/ctNBeFgOJxj1PosJ1c5rTxxxTHGZFmx2ayDnwoAqUePgD5+/PgVV1xRX19fV1e3aNGivLy87OzsK6+8srGxUYfSU8f7B2uIqKvMxogYERGZg8p9D65PXImRSMRkMmP8BkCa0uNX97bbbrv88suLioruvffeK6+8MhQKBYPBJUuW3HnnnTqUnjr2v/keEYXzTdwH60BnNSsXzKhIXImCEML6ogDpS4+96Wpqav71r38Zjcb33nvvt7/9rcViIaJvfOMbBQXja+vSA28cVE1c7whoIpohJbZ3WBSlvDzsnwKQrvRoQV966aXPPfccES1atOiVV16JHnz11VdLSkp0KD11HD1yIlQQG1DBaWzhxbMSV5yiKIxp2D8FIH3p0YJ+6qmnVq5c+f3vf7+wsPCGG25YsWKFpmnV1dV/+ctfdCg9RTDG/L5g55wsImJEHJGjXZnx4crElej3+6xWzE8BSGN6BHRxcfHBgwf//e9/Hz9+/PLLL8/MzCwrK1u+fHmiV0BOKbUn6oko4uDpgw5oZ0fk6o8uTVyJoijk5xcm7voAkGg6RSTHccuWLVu2bFnvEa/Xu2/fvmuuuUafCiTdn3+zg3Gkmj5YpJ9RYbuauAWMZFkm4sbVn0CAsSdpv8DV1dWrVq0SBOF8Jzz33HM///nP+x08duxYRUUChz0kzj///rpi43u38TYLapYjgY/vAgGfw4H+DYD0lrSAXrRo0QDpTESrV69evXp1v4P33XdfS0uC594lRsAflPqsAW0KqtMvnJa44iIRCf0bAOlO7ykMgUCAMaZzoUn337cPE1HYHftzaAmot92xJkHFqapKxPM8xm8ApDc9AloQhE2bNlVVVVmt1oyMDLPZXFlZuXHjxkgkokPpqWDbU38mIn9xbJF+e1CbcWFVgooLBoPRweYAkNb0COgNGzbs379/69atLS0tsiy3t7dv27aturr67rvv1qH0VPD+O0cYRxFnrElbqCWwg1gQwg6HM3HXBwB96NEHvWPHjvr6+t45x263e+HChfPnzy8vT+A6QamDMRYICqL7jFv9scsXJ65ETdPMZiwxCpD29GhBl5aW7ty5s9/B3bt35+bm6lB60h05fELT1GCfOYSWgPqJWz+WoOLC4bDZbErQxQFAT3q0oDdv3rxmzZpHHnlk+vTpLpcrGAzW1NR4PJ7t27frUHrS7X71TSIKZ8dCM6NdTtwi/eFw2G53JOjiAKAnPQJ6wYIFtbW1u3btqqur6+zsdLvd69evX7p06ThZJmL/m/8loog99s9KWSiBHdCSJObk5CTu+gCgG53GQRsMhhUrVuhTVqo5daKBiGRH7K/RZbMvTFBZgiAYjUaO4wY/FQBSHpZyTyxNY6FASLHF7rNB0pZdeWmCiguFgphACDBmIKAT69Sx00QUyjH1Ts4xCdqiy+YmqLhIJOJwYAdCgDECAZ1Y7793jIi8Ey29nQ42kazWhMwiEUWR5zn0bwCMGQjoxPrn3/5NRLIj1tdfFknUCnaiKDidGQm6OADoDwGdWMeOnCQitc+45LnZpQkqSxRF7EAIMJYgoBPL7wuqJk4189SzkzdbsnheIgpSVUVV1XEychFgnEBAJ1BnRzcRiVk9/RsckTnArv5IQnZRCQaxgTfAWIOATqB9e98lorA71sHhCKoZGQlZxigUCrlc6IAGGFMQ0AkUHcIRzosFdL6UkDWMVFXhec5kwhIcAGPKQAEty7Ju9RiTjh4+TkRiVqxfeHpOQpbgEEUJ/RsAY0//qd6RSGTr1q2vvPLKe++9197ebrPZLrzwwksuueSee+4pLMQWSsNTd7qZiFRTz19BRnT57FmJKCgUCmVmJnCHQwBIijNa0M8888w111wTDAbvv//+N954w+v1njx58nvf+9706dPXrl374IMP+v3+ZFU0Hfm9fiU2Q4UsQXXenBlxL4UxTZYjWAAaYOw5owXtcDheeeWVvkO17Hb7rFmzZs2adcstt7zzzjvt7e0ZGXgSNSSdnm4iUs2xm2kJqDNmxX+bK0mKWK1WTCAEGHvOaEGvXr3aYDC8+OKL/U569tlniejiiy+uqKjQr2ppbu/uA0RnbHPlEshsiX87VxDC6IAGGJPO8ZDw61//+g033NDU1ERE77333mWXXfaTn/xE94qlvT279hFRV1ls2Q2zNSG9EJIkWSzWRFwZAJLrHAG9f//+K6+8cvHixTfddNOHP/zhdevWvfnmm/rXLN0dOXySiKQsQ+86duWuvLiXoigyEcMEQoAx6RwBzfO82+1mjHk8nuzs7IqKCvRvjoDP6ycixWLovXerF8R/o9hwOGy1YgFogLHpHDuqLFu2TFXVl156aebMmbt37/7c5z53ySWX/PrXv9a9bmlMEiNCWOzbAW2QtGWXzY9/QZKYlZUd98sCQCo4Rwv6tttu271798yZM4no8ssvf/fdd/FscLg6O7qIyF8U6xo2SWzS5OL4lsKYJssyJhACjFVnBPRrr71GRGvXru3bp2E2mx9++GEiOnHiRF1dnc71S1N7/r2fiMSs3tvLMoX4z6oPhwWzOSFr/wNAKjgjNQ4dOnT11Vc//fTTR48ejc7zZow1Njb+/e9/X7NmzQ9+8AOnMyEL/Yw9/371DSJSYoOguSmB+A+0kCTR4cB3BGDMOqMP+gtf+MK6deuefPLJu+66q7q62mKxhEKhyZMnX3LJJd/5znemTJmSrFqmnebGNiLS+vQ9lOcVxLcIxlg4HHK70QENMGb1f0jodDofeuihhx56iIg6OzvdbjfPY8W7YQsGwpqJEzN7bi+vsOtWXhHfIqLDnzHABmAMGyh8c3JykM4joKma1+sPZ8faz7zKFiyaE99SBEGw2TDADmAsO8cwu1deeeXRRx/1er19Dx46dEivKqW907UNHJHcd5kkn5LpjvMaJrIsYQU7gLHtHAG9bt269evXf/zjH8f8tJFpbW4nolB+bGJ3gZeP783UNFVRVPx/AzC2nSOgBUH46le/is7NEas5fJKIZFssPSc5c+JbhCCIdjv6NwDGuHM0wVasWLFlyxbG2NmfgqF46a+vEJFRjt3Aypw473UQDAawQBLAmHeOgK6vr7/rrrtyc3OnTZt2wQf0r1n6am/tJKJAYU8XhzmoFWRlxbcIRZGxxCjAmHeOLo4tW7boX48xQ1EUSYr07d+IOPnlH4rnMkmiKKL5DDAenBHQubm5Ho9nxowZRPS73/3u5ptvTlKt0ti+Pe8SkZgVu7HmoDJ3wew4FhEKBS0WzPAGGPvO6OLw+Xy9H9922226V2YsqDlygogUa2zMRrYvzo9bJUnCDG+A8QDjtOLs1PF6IrL4ld4jM/l4DuFQFMVgMGCAHcB4gN/zOGtv7SCiwITYIOh5s+K5k7coing8CDBOnNEHzRg7ePDg2R8T0ezZ8exFHcPqT7cQUaAgFtCFFRPieP1wOJiREecxIQCQms4IaLfbvWLFirM/JqKOjg5d65WeNI11d3kZTxFXrA96cVXcBikyxhRFsVoxhANgXDgjoJHCo9TS3E5nziHkZVYxpSxe149EsIE3wDiCPuh4amlsISKhzxg7R5iZTOcYbD4yoihareiABhgvENDx9PZbh4go1KcDOlc1n//0YQuHwxgBDTB+IKDjqb62iYiozyomFZa47XiiqqrRaDAa49YeB4AUh4COp+NHTxORbI/d1WtmXBSviwtC2GjEBt4A4wgCOp66PF1EpBliUweLS4ridXFRxBaxAOMLAjpuZFnx+0NEpJp7AtoSVOfPnRmXizPGJEkym9GCBhhHENBx4+v2EREzchFnTzex5DTkFeTG5eKiKJrNZiLsogAwjiCg4+ZYzWkikhy9t5Rl+Inn4xOpkiQ4HI64XAoA0gUCOm6OvH+czlholMticeuRCIcFjIAGGG+SE9BvvPFGUspNqKaGFiIKTIiNU55gj8+u27Is8zyPFewAxpvk/M5fffXVSSk3oepONxGxiCO2CkdVYXFcrhwKhbD+BsA4pEdAf/GLX7SdKRgMRj/QoXTddLR5iDj6YLNdRnTt7HlxubIoCk6nKy6XAoA0ose0tE2bNgmCcPjw4c2bN+fn5xNReXn5qVOndChaT63NHUQkZfTcUo5o9sxpo7+sqqqMaZhACDAO6dGCttvtW7Zs+fKXv3zzzTcfOHCgoKCA5/mCgoKCggIdSteHpml05jp2ppCamRmHZm84HLZYxtS/GgAwRPq1y6699tr58+evXbt2x44duhWqm6bGViLyToz1FBs53u6IQ7BGIhL6NwDGJ10fEhYWFv7tb3+78MILV61apWe5Ojhy6DgRMS62TlJuyHD+04eOCYKAFewAxie9R3FwHLdhw4Znn33W6/Xu3LlT59IT58TRWiKSnbH/SCaH4zDuQhCiEwgBYDxK2qOn6urqVatWCYJwvhPee++9N998s9/BQ4cOpebYj2M1tUTEq7EWdHlxHJZJEkVMIAQYv5IW0IsWLRognYnIarW63e6zD6bmfI221g4ikq093RrmoLJgzvTRX1aW5cxMbBELME7pHdCBQMDpdHLc4CtUVFZWVlZW9ju4d+/elpaWxFRtVMIhkXEULOiZ2x1xGmdfPNqA1jSNMZaaf5AAQAd6/PILgrBp06aqqiqr1ZqRkWE2mysrKzdu3BiJRHQoXR8tTW2yPfZU0BBhefmj3UslHA6aTFhfFGD80iOgN2zYsH///q1bt7a0tMiy3N7evm3bturq6rvvvluH0nXQ0d5JRJoptteVq1srLMof5WVFUbTb0QENMH7p0cWxY8eO+vr63od7brd74cKF8+fPLy8v16F0Hbx74H0iklzG3o4boz8ylG6cgUmSlJuLAXYA45ceLejS0tKzR9Tt3r07Nzc+i9knXX1tIxEZJK33SKlqH+U1JUkymy1YoR9gPNOjBb158+Y1a9Y88sgj06dPd7lcwWCwpqbG4/Fs375dh9J10N7aSUSBCbEByyWTR7uOXTAYsFrRfAYY1/QI6AULFtTW1u7atauurq6zs9Ptdq9fv37p0qUGQ1zm2iVfZ0c3Ecl9Fhq9ZMYFo7xmJBJxu0f7mBEA0ppOw+wMBsOKFSv0KUt/zU1tRBTps1LSvAtGtY6dqqpYoR8AEAFx0N7mYRxFXLG/dtmuUS1vJAhhrL8BAAjo0QqHhIA/1HehUXNInVJRNpprhkLh1JzRDgB6QkCPVnQl6EhGrPlsDmo2+8jbv4wxRYmYzWhBA4x3COjROnrkFBFpfGyZpInKqMbYiaJgtdpGP4waANIdAnq0PO2dRCRlxcbYuR2j6oAOhwWbbbTDqAFgDEBAj1ZjQwsRKZZYg7fcOar15xRFxh7eAEAI6NE78OZ7RCS6YoOgC/PyRnw1WZaJCAPsAIAQ0KPX0dZJRL1L2TGi8oLCEV9NFDF+AwB6IKBHy+cNqGZOsfbcSY5o7tyZI74aBtgBQC8E9Kgoiubz+iN9JnlbAmpRyQhb0KqqMqaZTNiEEACIENCj1NLUSkQUW8aOsurFLHfGyK4mSaLFgseDANADAT0q0VU4fJNik0oyM0eYzkQUCoWwQj8A9EJAj0pLUzsRSX2mEWZaR96DrGkaluAAgF4I6FHx+QJ05jp2lxZPGdmlIpEIz3OYQAgAvRDQo9JU10JnrgQ9a6T7eEUiEasVEwgBIAYBPSqHDh5RLLF7aBS1zMwRzvMOhULo3wCAvhDQoxLwh/q+tHuVvIIRboPCmGo2Y4AdAMQgoEelvdXDeCKi6Fp2GfViyaSR7EYYiUhGoymeNQOA9IeAHrlo8zlYYKYPNt/W8u0TivJHcKlgMIQFkgCgHwT0yAX8QSISM2NPCNVAZGSXkiQRS4wCQD8I6JHz+wL9jlS4c0dwHVVVeJ4bM3ucA0C8IKBH7q03/ktEzjYl+tLWrVTaR7LQqCAI1lFMbwGAsQoBPXKnTtQTkbesZ2yc4DaaRrSOczgcttnQAQ0A/SGgR669uYPOXKq/fOawpxEyxhhjJhNGQANAfwjokevu8mkGrncaoUnQZs+oHO5FIpGI0WjCDG8AOBsCeuQiEblvrJpDSn7BsB8SCkIYA+wA4JwQ0CPEGGuobw67Y+vYmcJscsXE4V5HkiQENACcEwJ6hHzdfiLiWeyIUdRMJuN533AuqqoqiowBdgBwTgjoEWpuaicig9SzmYohwiYHLMPtSg6FgpifAgDng4AeoV07/0NE3WU9vROqmZNyhr3UUTgcdjpHuPodAIx5COgRCgRDRBTKjfVp5AxztypN0zRNwwp2AHA+COgRam1uJyJDn7U3Vlwyb1hXEEUBjwcBYAAI6BEK+EJEpBl7Op3NQXVOZdWwriAIAjqgAWAACOgRamlqJY6E7J4ujojTkD9hGIOgGWOiKKIFDQADQECPhKZqHk83MbJ6e1ZKymiSJk4qGvoVIhHJYjFjAiEADAABPRKiJBGRbOPFrJ4WtFHQrNZhrKcRDoexgh0ADAwBPRKnjtcRkWyLTTAxZNl4wzBuZjgcttvRAQ0AA0FAj0Snx0tEoZw+Y+z4YTSfFUU2mYw8jwmEADAQBPRINNY1E1GwKBbKLscwBkGHwyGzGY8HAWAQCOiROHm8johUS+wR37zcYSyTJIqi0zm8WS0AMA4hoEfC7wsSUcQe66OYO7F8iO9ljKmqZjSaElIzABhDENAjIQpi33Q2RLRJZSVDfK8ghM1m7J8CAINDQI+E1+vvnaJCRCZBmzF7qNMIw2EBOxACwFAgoEfi2JFTzBDrgLa3y6FgeIjvlSTM8AaAIUFAD5soSkQUccZunatdKijMG8p7I5GI2TzsZaMBYHxCQA9bdB27YH6ffmTG8fyQMjccDlss6IAGgCFBQA+bp72LiMTM2EPCiqwhNZ+JKBwOORzOhFQLAMYc/QLa6/X2O+Lz+XQrPY4a61s0U+y+8TKbVFgwlDfKssxxhB0IAWCI9AjoQ4cOzZgxIycnZ8qUKdu3b48elCQpKytLh9Ljru50IzGt96XNpxRMyB/KGwUhbLNhfgoADJUeAX3nnXeuXbtWFMVf/epXd95553/+8x8dCk0cvy8Y7rP9oDGoTiwb0kKjgiA4hjMjHADGOePgp4za4cOHX3/9dZ7nly1btmXLlvXr1x88eFCHchOkqaGV9XkiaBZZSengAc0YY4yZTJhACABDpUcLetKkSa+++mr042uvvXbWrFn33XefDuUmiCwrYlYsZ80heSgtaEEQDAY9/hwCwJihR0B/5zvfuf766y+77LKuri4i+vnPf75///7FixfrUHQidHf5eDXWB20Umck8eLs4GAxggSQAGBY92nQrV66sqanZs2dPdAhwVlbWnj17tm/ffuDAAR1Kjy9NY411zfyknrna5qBSIJpLJg7SgmaMybJss2ELFQAYBp3+6S4pKbnpppt6X5pMphUrVrhcrgHeoihKIBDod1CSpITUb8i8XT4i8pf0PCSMOI3din/QWSqCIJjNZiJMIASAYUhar2h1dfWqVasEQTjfCS+88MLmzZv7HTx27NjUqVMTXLWBtLa0E5Fi7RnLzIjyMwcfLCiKIacT81MAYHiSFtCLFi0aIJ2J6OMf//jHP/7xfgfvu+++lpaWRNZrEB3tXYwjyRWbbFKUP/g0wkhEycrCCnYAMDx6T/UOBAKMMZ0LjaMuTzfjqfcLsHmVkokTBn5LJCJxHPE8ZtUDwPDokRqCIGzatKmqqspqtWZkZJjN5srKyo0bN0YiER1Kj6/GhlYx09jbl2wUtQklg8zzDoXCdjv6NwBg2PQI6A0bNuzfv3/r1q0tLS2yLLe3t2/btq26uvruu+/WofT44jlOyoz1CxkirKhokHne4XAI4zcAYAT06IPesWNHfX19b0i53e6FCxfOnz+/vHyo+/iljsb6ZrXPSkk2rzJj1kB7qchyhON4oxFTVABg2PRoQZeWlu7cubPfwd27d+fm5upQenx1dnp7nxAyIocn0tXZPcD5oVAY81MAYGT0aNlt3rx5zZo1jzzyyPTp010uVzAYrKmp8Xg8vSvbpZFOT3f4gp5B0BwR42jylIkDnB8OB/PzC3WpGgCMNXoE9IIFC2pra3ft2lVXV9fZ2el2u9evX7906dK0WxlZ01jD6Wb5oljD3+JXs3PPOw5a01SOM6B/AwBGRqfsMBgMK1as0KesBGKMiCwBNdrL4WqNcIwGGDQYCPgtFvN5Pw0AMCAMzh2GUEiIOAx9Z6kUlxaaTOf9IxcKhZ3OgaazAwAMAAE9DCePn6Y+y25wjHJy3ec7WVUVTdPMZrSgAWCEENDDYDKbwu5Ye9kYVC2W8y406vf77XaM3wCAkUNAD8PJo6eNUp/dCP1KQdF5pxFKkjjwcn0AAANDQA9DY0NLODvWgrZ3KxkZ524jK4qiadjgCgBGBQE9DPWnm/0lsUXpIjY+r+Dcc23C4TD6NwBglBDQw3C0+kTEEbtjFp+cl59zzjPD4bDdbterXgAwNiGgh6GpodUc7umDtgRUY4SKis/RB61pKmMYvwEAo4WAHgbFwkfsPXfMIGlEVDLpHItBh8OC1Yp0BoDRQkAPlSwrxohmFD5oQYfU7OysLHfG2WcKAjqgASAOENBD1VjXHM42KbYPWtAi6+rynn0aYywSiZjN2OAKAEYLAT1UsqJoxj7bcvOsvGLS2adJkmg2mzkOG3gDwGghoIfq2JFTap9hzc5WOTc/++zTAoGAzYbxGwAQBwjooWprbhezYgmt2PjCCefYz1tRVIcDAQ0AcYCAHqrWZo+9S+59affIE4r6B7QsR3ie4zjcVQCIA0TJUB2rOekvig2e42WWf1YL2ufzY39YAIgXBPRQKYqiWGMrQWsmrri0/15WghB2Op361gsAxiwE9FD5fEHV0HO7zAHVKGoFhWe0oCVJMhqNPJ9m+3gBQMpCQA+JoihtbZ2iuyd8Iy6DK8NZVHLGPO9QKOhwYH4KAMQNAnpIAv4Q44mXe/YftHcqAX+w32BnQQihfwMA4ggBPSTBQEh0GzVTTyIbJa30zFU4JEkyGs3o3wCAOEJAD8mJY6dNIbX3pbVbzsk7Y5aK3+/D+A0AiC8E9JB0tHeF8kxEPV0cppCalXXGMkmKImMDbwCILwT0kLQ0talmnqini0Mz88UTY2PsJEnkOAPP42YCQDwhU4bk5PE6m1fmFUZEprCa2SQVFMTG2Pl8PozfAIC4Q0APScAf6iqzRVezk+0G1cwXfTBL/F8f8wAACgZJREFUhTEmiiLGbwBA3CGgh6S1qU12xkZoMEa5ee7ox6Io2Gw2rC8KAHGHgB6cqqqdnm4hy/jBAWYOq2VTJkZfhMMh7J8CAImAgB4KjogsfiX6wtatEpHV0rNwUjgsYANvAEgEBPTgmhtbVTMnZfS0oI2iNqEonzfwRCRJoslkQv8GACQCAnpwrc0dsi12o5iBcgtzoh+Hw4LLheHPAJAQCOjBNTe19zaficgoMKOh54FhKBS0WjGBEAASAgE9uJamVj42zZtsXnlyeSkRSZJkMBgMBqy/AQAJgYAe3OmT9cRY70tXSySvIIeIAgEfhj8DQOIgoAfn8wYDE2KbXUlOQ+mkIiISBNHhQEADQKIgoAfX2tJm9Sq9Ly1BtXRSkSgKJpMJ628AQOIYBz9l3GtqaHM6DdZuRTPzuUdDRlHLK8jx+XyYnwIACYUG4ODc2VmtMx2i2xhx8M0Xu4jIleGQZRkD7AAgoRDQg+ju8nV3eZU+46BzJ+SqqmKxWDA/BQASCgE9CF+3n/EkZvb0BZkDamlxQTAYwONBAEg0BPQgJCmimmJ3ySRqiqJKUgQbXAFAoiGgB1H93jGDpNm65ejLzAYpOzfLarWifwMAEg0BPYhjx2rFLKPgNkVfakbOaDS4XBkDvwsAYPQQ0IMI+EJ9V0qS7XxeYbbVak1ilQBgnEBAD6KlqZWZY70Zzg65atqUJNYHAMYPBPQgNFUVMk29L1UjlVeWJa86ADCO6B3QgUCA9Vl4KPUpqiZbYy1oe5eSmYEOaADQgx4BLQjCpk2bqqqqrFZrRkaG2WyurKzcuHFjJBLRofRROlp9UjXH7pLL6ZhQnJ/E+gDA+KFHQG/YsGH//v1bt25taWmRZbm9vX3btm3V1dV33323DqWPRqenm4hUS89dMgma0BHAADsA0IceiyXt2LGjvr6+d2aH2+1euHDh/Pnzy8vLdSh9NHieZzwnuXqW5JdtfHnFpORWCQDGDz1a0KWlpTt37ux3cPfu3bm5uTqUPhrHjpzkNGYOadGXrpaI2WIa+C0AAPGiRwt68+bNa9aseeSRR6ZPn+5yuYLBYE1Njcfj2b59uw6lj4YsKxGnIeLo+TNmFpQJRXnJrRIAjB96BPSCBQtqa2t37dpVV1fX2dnpdrvXr1+/dOnS1N/Nz+cNRJyxSqpGw8TJJUmsDwCMKzot2G8wGFasWNH3iNfr3bdv3zXXXKNPBUbm+NFa0RULaFNYdTjsSawPAIwrXLJGJe/du3f58uWCIJzvhBdeeOE3v/lNv4MHDx7c+I3/4w08EdUea25u6MgrdFfOmBj9bO2x5saG9hNF2hXlFdEje+prT3V6ijOzlpdPHcGRKydPjWiqSMrfmk8c83uWezLXfmR5tBu6vLx82rRptbW1R44cib4LR3AER8bhkQMHDlx88cVTpiRkgnHSAnpQsiwHg8F+B//yl79IonzdddcSkdVqNRgMqqqKohj9bPSITwj15r7ZYomeE5GkERxpbW7/6pe/1zrV0pTLsk6Giw6Hntv5lCvDSUQOhyN6TigUir4LR3AER8bhkfvvv/+ee+4ZIwEdCAScTueIhxL/6U9/8ng8eg6grj1R//SWP3W0dxaXFq6766YJxQW6FQ0Aqe/2229/8MEHExTQevRBC4LwxBNPbNu2ra6uTpIko9E4efLkW2+99Stf+YrZbNahAqMxuWLixsf/J9m1AIDxCDMJAQBSFGYSAgCkKMwkBABIUZhJCACQojCTEAAgRSVtJiEAAAwMW14BAKQoBDQAQIpCQAMApCid+qDjJTMzc9OmTc8//7zO5e7bt89kwlL98ccYkyTJarUmuyJjkCzLRISf20QQRXHOnDnRuR0NDQ29kzziLnUXS0opV1xxxa5du5JdizHo+PHjjz/++C9+8YtkV2QMevrpp3mev+2225JdkTEooetv9IUuDgCAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSVZsPsksVoxI1KCJ7neR6thIQwGAy4twmi288thtkNSSAQcLlcya7F2IR7myCRSISIUn/TonSk2w8tAhoAIEXhPyAAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgB7EW2+9NXfu3Ly8vHXr1gmCkOzqpKulS5dyH/jYxz4WPXjOe4sbPnSPPfZYc3Nz78uh30/c5EH1u7dJ+wFmcH6KopSWlv7lL38Jh8PXXXfdt771rWTXKF2VlJScPHkyGAwGg0FBENh57i1u+BCFQqHnnnvOYDDU1NREjwz9fuImD+zse8uS9wOMgB7Iyy+/PHPmzOjHe/bsqaysTG590pQgCE6ns9/Bc95b3PAhuvXWW2fOnGk0GntDZOj3Ezd5YGff2yT+AKOLYyB1dXUzZsyIfjx9+vS6urrk1idNnT592mg0zps3z+12X3311adOnaLz3Fvc8CF69tlnDx06lJub23tk6PcTN3lgZ9/bJP4AI6AH0tnZ2bskisvlkiQpEAgkt0rpyO/3L1q06I9//GN7e/tFF130yU9+ks5zb3HDR2zo9xM3ebiS+AOMVTQH4na7g8Fg9ONAIGA0Gp1OZ3KrlI7mz5//0ksvRT9+9NFHnU5nV1fXOe8tbviIDf1+4iYPVxJ/gNGCHkh5efnRo0ejH9fU1EyePJnjuORWKR299dZbe/bsiX5sMBgMBoPFYjnnvcUNH7Gh30/c5OFK4g8wAnogy5Yta29v/9e//qWq6pNPPnnrrbcmu0ZpKRgMrl69uqamRtO0TZs2LV++3OFwnPPe4oaP2NDvJ27ycCXzB3iUDxnHvAMHDsyZM2fixIlr166VJCnZ1UlXP/rRjyoqKvLz81evXt3a2ho9eM57ixs+dIWFhX2Hgg39fuImD6rfvU3WDzAW7AcASFHo4gAASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoAIAUhYAGAEhRCGgAgBSFgAYASFEIaACAFIWABgBIUQhoGF+efvrpK664gog0TZs7d+5LL72U7BoBnBfHGEt2HQD0wxj7/+3Yoa2DABSG0ZC8EXDoupoKHGEBhsB3FDZA1KCZAdHCAg2SCTB41JuAJx834ZwJfvXl5pZl+Xw+t20bhqHv+7MXwSGB5nLmea6qKkmSaZqyLDt7Dhz6OXsA/Lf7/X673dI0VWeC84Pmct7v97qu4zh+v9+zt8BfvDi4ln3fH49H27bLsrxer8/nc/YiOOSC5lqapsnzvCiKuq73fe+67uxFcMgFDRCUCxogKIEGCEqgAYISaICgBBogKIEGCEqgAYISaICgBBogKIEGCEqgAYISaICgBBogKIEGCEqgAYISaICgBBogKIEGCOoXWX8tPr4ED3QAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "plot.ecdf(NONLEAD_tabf0502$Dist_SNP_to_peak_center[!duplicated(NONLEAD_tabf0502$Peak)],col=enh_cols[1],cex=1.2)\n",
    "plot.ecdf(LOCAL_tabf_2$Dist_SNP_to_peak_center[!duplicated(LOCAL_tabf_2$Peak)],col=enh_cols[2],add=T)\n",
    "plot.ecdf(CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center[!duplicated(CMLEADS_SNPinpeak_2$Peak)],col=enh_cols[3],add=T,cex=0.6)\n",
    "\n",
    "## the no-effect peaks are much further away from the peak center --> pre-sample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a6a1fa",
   "metadata": {},
   "source": [
    "### Pre-Sampling the no-effect peaks to be similar in SNP-to-peak-center ditance to the larger set of LOCAL caQTL peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "acfe8220",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 7625\n",
      "[1] 7625   25\n",
      "[1] 7611\n",
      "\n",
      "\tWilcoxon rank sum test with continuity correction\n",
      "\n",
      "data:  (LOCAL_tabf_2$Dist_SNP_to_peak_center) and (NONLEAD_tabff$Dist_SNP_to_peak_center)\n",
      "W = 27480507, p-value = 4.958e-09\n",
      "alternative hypothesis: true location shift is not equal to 0\n",
      "\n",
      "\n",
      "\tWilcoxon rank sum test with continuity correction\n",
      "\n",
      "data:  (LOCAL_tabf_2$Dist_SNP_to_peak_center) and (CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center)\n",
      "W = 7704154, p-value = 0.3046\n",
      "alternative hypothesis: true location shift is not equal to 0\n",
      "\n",
      "\n",
      "\tWilcoxon rank sum test with continuity correction\n",
      "\n",
      "data:  (CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center) and (NONLEAD_tabff$Dist_SNP_to_peak_center)\n",
      "W = 7518149, p-value = 0.007299\n",
      "alternative hypothesis: true location shift is not equal to 0\n",
      "\n",
      "[1] 7625\n",
      "[1] 7625   25\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3deUBU1eIH8HPnzsosMOyIIqAsAuIKKuKumS2mlWVZmfls0V+LLS8rnz19Zb323fz109KysrKeWpkv08RdNAMEWcVBBIVhX2af+/tjbMQ7w+rMne37+evMOZe551J9u5x7zrkUwzAEAADcD8/VHQAAAPsQ0AAAbgoBDQDgphDQAABuCgENAOCmENAAAG4KAQ0A4KYQ0AAAbgoBDQDgphDQAABuCgENAOCmENAAAG4KAQ0A4KYQ0AAAbgoBDQDgphDQAABuCgENAOCmENAAAG4KAQ0A4KYQ0AAAbgoBDQDgphDQAABuCgENAOCmENDgW/Lz8ymKGjJkiOVjdnZ2fHy8QCBYv369pWbr1q08Hq+kpKSzb3jllVcCAgIaGhq46C74NgQ0+LS33367pKRk/PjxKSkphBCTyfTiiy/OmjUrLi6usx958MEHdTrd66+/zmE3wUchoMGnWW6En3jiifHjxxNCNm/eXFRU9Pjjj3fxI0FBQQsWLHjvvfcuXbrEUS/BVyGgwfOcO3duzpw5ISEhYWFh999/f21traW+sbFxyZIl0dHRAQEBM2bMOHXqlKW+vr7+zjvvDAoKGjFixJEjR6zfM2LEiF9++YUQMnfu3NWrVxNCXnrppcGDB8+YMYMQsn79eoqiQkNDW1tba2pqZDIZRVE//PADIWTZsmVtbW3vvvsuxxcOPocB8ChNTU3h4eECgWDu3LmZmZmEkJSUFL1ebzabx44dSwhJTU3NzMykKCo4OPjSpUtms3n06NGEkOjo6FGjRvF4PEJIYmIiwzC//PLLsGHDCCFPPfXU6dOnCwsLCSHz58+3nMhsNk+YMIEQ8tJLLy1fvpwQMmfOHEuT0WgUCoUjRoxw2W8BfAMCGjzMm2++SQh59NFHLR8zMjIIITt27Pj1118JIXFxcUajkWGYu+++m8/nb9iw4ccffySEDB06VKPRMAxz3333WQOaYZjrr7+eEPKf//yHYZivv/6aEPKPf/zDeq4zZ84IhUJ/f3+xWKxQKCorK61NCQkJNE2bzWbOLhx8EIY4wMPk5eURQt5//32KoiiKOnz4MCGkoKAgPz+fEJKZmUnTNCFky5YtBoPhgQceOHnyJCHkjjvuEIvFhBBLQNtVV1dHCImOjrbWJCYmPv/8801NTVqt9pVXXomMjLQ2xcbGmkymxsZG51wlACGE8F3dAYDe4fP5hJAHH3xw6tSp1sqEhITffvuNEBIUFMQ6vqWlpeNHPz+/zr5ZKpUSQurr6ztW6vV6VsHCkuZdfBvAtcMdNHiY5ORkQohAILjzzjvvvPPOpqamEydOUBSVlJRECPn111+NRiMhZOHChXK5/IsvvkhISCCEbNu2TafTEUK2bt3a2TcPHDiQEFJaWmqtycnJee2114KCgmQy2cqVK8vLy61NpaWl4eHhIpHIWdcJQPCQEDxNU1NTSEgIj8dbsGDB7NmzCSGDBg1qbW01mUzDhw8nhAwbNmzy5MkURUVERKjV6vr6esutcVxc3Lhx4yiKIp2MQWu1WolEMnnyZEuT0WgcNWoUIeTtt99euXIlIWTGjBmWJsu8EevjRAAnQUCD5ykuLp41a5ZSqQwNDb3vvvsuXLhgqa+pqbn33nsjIyP9/f1nzZqVn59vqT9w4MC0adP8/f1TUlI++uijzgKaYZhHHnmEpmmVSsUwzGuvvUYIiYqK0mq1TU1NlsGTTz/9lGGYd955hxCyb98+bq8bfA7FMIxL7twB3FBVVdXgwYMfffTRf//7350dYzab4+PjY2JiLPNGAJwHY9AAV/Tr12/ZsmWffPJJe3t7Z8f89NNPZWVlL730EpcdA9+EWRwAV3nuuecSExMbGho6m6GhVCq//fbbMWPGcNwx8EEY4gAAcFMY4gAAcFMIaAAAN4WABgBwUwhoAAA3hYAGAHBTCGgAADeFgAYAcFMIaAAAN4WABgBwUwhoAAA3hYAGAHBTCGgAADeFgAYAcFMIaAAAN4WABgBwUwhoAAA3hYAGAHBTCGgAADeFgAYAcFMe9tJYtVq9b98+V/cCAOAyHo83e/ZsgUDgjC/3sIDeu3fv9u3bJ02a5OqOAAAQQsjmzZuHDx8+aNAgZ3y5hwU0IWT8+PEPPvigq3sBAEAIIcePH3fel2MMGgDATSGgAQDcFNcB3dLSwjAMxycFAPBEXAS0RqNZu3ZtQkKCWCxWKBRCoTA+Pn7NmjV6vZ6DswMAeCguAnrZsmXZ2dkbN26srq42GAw1NTWbNm0qKChYunQpB2cHAPBQXMzi2LlzZ0VFhUQisXxUKpXjxo1LT0+PjY3l4OwAAB6KizvoAQMG7N69m1WZlZUVHBzMwdkBADwUF3fQ69atmzdv3qpVq5KSkuRyeWtra2FhoVqt3r59OwdnB/B2jMFg7GwlG8MwRmOnrd5EY9DXtbf2Uyh5FOXUEzEMQzn5FFZcBPSYMWPKy8v37dunUqnq6uqUSuWSJUsmTZpE0zQHZwdwLct/zwxj1usNFEXpdDqKYnQ6A03zaJo2mxmdTsMwRCLxk0ql9fV1DQ0NZrOJoiihUERRlF6vMRrNhBAejxaLRe3t7YQQoVCk1+sIIXy+wGg0WE4UFhYeFHTVX6U1NZfU6lpLecCAKLlcQQgxM4zGoFc11PqL/UpqLzZp29OiBoXLAlTl5+WBip/LcssuVJ/PPdti0klo/jD/yAC+JHJA2Oixw4JDAgvyig/9nq1WN1CEUvhL5Qp57SV1WYkqNDxo6szMCVPGEELUtfXbvvq59lJd0tC4m2+7TiC4EjJZe499/s2OkvpLdHzQ+PGjnpl6s1x0eeTzt90Hf/vlIEVRmZPTzpVV5ucVBQYp775/TmLy4LIS1acfb61TN/SPili89C6NiDz/01c7808SQl6+Yf7jE2cxDPPiL9++tf8ny1d9v+ip6xJSu/gncrZU9eqqD8+WqZqbWjMmjn7mH49EDgi3e+TJ43mffbz12OFTGRNHT7lu/Neb/lNWoiKEDE6I+fuqR4aPSu7dvwq9R7lq0ltjY+OxY8dmzpzZq5/65ptv1Go1ni4CZxjGrNFoKIpiGKLXa7VanVarIYSIRGK9XmcymQUCvslk1um0ZrOZogiPRwsEQoGAbzabWlvbOO6tRk80Te0//PBfrcE4OHHgzGlpHVufeuiVOnWjmUcRiuEZL1e2Bwlog1nUbLJ8NIl5tNZMCEMI+yZx0rSx+3872sXZ73ng1jvvnX3zlPutNZmT09/4aBWPRxFCtnz6/bv/3mBtqh4hH3f92K33PU5R1Bcbv3/vtQ22X0gIeWPdi08/str6MXXEENVY/yPVpdaa929dJKDph7/9v44/VfrCu+HyALtfWFtTf+PEezvWDEmJW7f5FT8/CevIc2Xn77jx4S6u9+cDnweHBP7tb3977rnnnLTU22UBffjw4WnTpmk0ms4OOH78uO2+SEeOHJk/f354eDghJCYmZuDAgRcvXiwsLLS0ogY1fa6hKCoyMlIiETU2NjY3t1iO8fOTCAQCo9HY1tbuETVNjS25JwtrLzUqg+RpGUNlcmnHYypUVe1t7YSQ6kp1Q12zwl/aPzrM8j2OqhHwhTkni4rzz46ZdPkedtiIlMlTJ1RXV3/z5X8sNccP5JYUVkSMiLrp1kn9AwL5NP3Ze9+VFFZERYdPmDG64zFR0eEzZo/XanQda6bdMZGWCAghP58rzKo6O3tw6pR+0Tqj0VozOnTA/clpApomhFwy6usZYwCPH0ELLd/8Z/7ZXzfvYZ1LPjwqLj6qv8TPUnNRo1Gb9K0FVdHBSlZ/Ov4UifBPHBorp8w3XzfD2/biyMjI6CKdCSERERGjRo1iVR4/fryysjIpKYkQolAoCCEymSwwMNDSihrU9LCGYcwmk6G4uNBoNMrlMssxGk27Xq9lGCIUXv6PmcejCSE8Hs+mhnbPGq1G36Bubm/TEEJ0WkNgkLDjMeqL9TWX6gghFypqmptapTKJ0Xj5RtpRNVVVdY3qpvr6pnOlFyw1e7P+WPnv9QFyv7iQUEtNY0MLIaT9XMOZE8XlAiFFXa5pamqz/pS1prjgnNlk7lhTUnCOSPiEkOazVUF1mqrG8qI4ojUYrDWtly6WCy+I+QJCiLquqbVd1yoWGEIuR231mUrbcwnK65toiUB5+d+E+qY2nUbfWFF3rrGd1Z+OPxUVGSKlxBerzhGncdkddN8sX768urr666+/dnVHwJO0tbVdunRRp9MyDENRpPN/5e38Xe9ZVj7xVmXFJUs5MTl2xb8e6ti67L7Vba3tTu3Awgfv2PS/33Ss2fz9e4lJgwgh6Yk3dqxvCxGeG69Q/+sTiUD4yMLnTh7LtfuF06+fsOeXAx1ryicGaAKv3Fn+a9Ydqf0G3rLh9Y7HVK/+2Dq6zfJHdt7D965gVe7c91lYRAir8uSx3EcWPmf3Syx27PssPCLEqUMc2IsDvFZTU9OZM2cKCk6rVOVarcZyL9LlDYm7p3NL8+VB7eamVkuhtqbB2vTJ+1vP/5XOhJDC/LMfvbmlMP8sISTn5JmVy9+5xnQOt4kwlkefXnT/g/NGpKVYa+68d7YlnQkhH376cseDa4b4/bDoKYlASAhZ9uRCu184flLaqleWz7njemvNQ4/d88uql8bHJFg+Pjhu+mMTZk2LS9lw58OZsYmEkJuSRp54cm1n6UwIGZk2dMn/3NOx5uW3nrVNZ0LIqDGpj/19sfXjnHlXutGvf9i7n6zp9hdy7bi4g66qqmppabHblJCQ0Kuvwh00dKumpqahodZkcsc/DRkzQ/Eos9lsMplMRnO9ukmv1wtFAr3OmPtH4dEDf06cljZiTLJEIm6oa/rtl0MXKi7J5H4SidhP5nfy6OmmxhapzI9hzO1tWrFYpNcbzGYzRVE8Hs9kMvW5V8r+QQ2VdZZyeEJkS3NbW3UjLeDRQoGIFgiFgsgB4bfMmznt+gnff/1z1t6jbc3tUrkkLCK0X2SYXm+srKiMjIq47oZJSUPjCSFGo3H/b0cvVatThiWkjhjS8UR16oajB05ebGuOSo0ZlzBEJhJbm1qa204cyyGEpI0d1tbaXlhQFhSsTE6Nt0xoO19RdaHi4qC4gSFhQZbja1ubFWI/Eb+Pg7QNdY2nc4rEEtHQEUPEYlEXR+q0+uqqmn79w4RCgdnMqGvqhCJhgFJhPcDjHxK++uqrK1eupCjKdmVKdXV1r74KAQ2d0Wg0lZUVBoPBSd+vN5n0JqOApilC2nX6/LxS2shIRMKCvLNmkzG8f5ifn1hdW597soimeZUVF7UaXUh4kF6n17Zr29o0hFAURZkZhqO7dIoQAc1XiPkSQXCQUtBubm1uoyWCAcNiEhJib7tpRmttI18oDAoKsGSN2WQ2mkxCofdPl3Y4pwY0Fw8JV6xYMXjw4NWrV+fl5XFwOvApZrNZra6tq6u1d6fRxzFliqJ4PLpe15alKj56ofy0urqxoiE8p0XcaqLMDOluCNt61urKGutXEsuE6O4713mfKcoopGjC8+PzFQpZv4jQzKljR6QlF7bXnqq7IJWI7xqeQShqQEAgn9eDFQYhV90t8WiekMaAp9vhaBbHrFmzsG4QHMtkMpWWlphMxs4P6T6deTyKzxfy+Xx/f0VAQOAxVenXfx7+tSBXc6Qi8JyGZ2QIQ5SEKHvUI8ruWa2Jy4pe284FKBV6vUGvN/CFgrABoSGBAUGBAdNnTRg9ZhgtFhhNpo5jAlYpJPH2HnUPPA9HAS2VSj///HNuzgVer729taLivNnc61FXHo8IhWK53D8gwF8gEBJCjGbTkeKid9/bULHvNDFa7o4ZP0JJbX62zzM8KJvCX52hRGKRTCGdecOkZU/eT/O7u+3lY/zB53jeOwnBl2m1mnPnzvUmmhlCKLFYPHBgjHVrAZPZfOBs4ftvfnp+/xliZg9YMMT+Pgt9HjvmUZRQLBIIBUHBAWljhg8dmThydHJouNMnAIAXQECDB2AYRqVStbe39vxHeDxKJlOEhISKRFee0f+0//Dqp18nLV29KeIaH+LxKF5gSMDceTPTM0cMiouWyWzvxQF6CgEN7u7s2TLL9hc9JBAIBg6Msa61I4TU1zffPWdpfW1D58/3+oAhhBJJhEFBgdGD+s+bf9P4KWnd/xBAbyCgwU0ZDAaVqlSv7+loBkXx/P39IyL6ddwK8stvfn53zXrG2MWDxF6gCFEEyMdPSUsfM3zarMyO9+YAzoCABrdjMOjPni3r+coLgYAfGxvXcffaIwf/ePGZNxsbGi0fOzzf682jPoqiebyo6Mjho5NGjR02/bpMHiaiAbcQ0OBe6upqLl2q6f44QgghIpEkNjbWestsNpsfvm9Fzsl81pxoym7RFo/nHxHw0EN3T58xXib34/d1lRqAo+BfQXAXDGMuLi7q4Y2zSCSKjo613DUzDPPyqvd/2varyWwmvZwPR1G8iNSo9R+uCQsO6mO/AZwGAQ1uoaqqsrGxsdvDaJrq3z9aKr0yNeKpR9Yc+P0YYZjO1onYxYh4Dy1fuHjhbZy9uwigDxDQ4GIdX8vUBR6PFxMTK/prKZ3ZZH5lzQc7tu7+azCjRzlrFvJSH7vulTsXhskV3R8N4GoIaHAZs9lcXFxoNpu7Pkwg4EdHx1oW/hFCfvz+v6/+80O9vtuJGVcNdWgD+C9/8NzNo8deU48BuIWABteoqrrQ2NjQ9TEURScmJlDU5bkTe/978IXlr/V4dgdFCNUeIjifLmOE/PV3LL55BNIZPAwCGrjW1taqUp3r9jCZTBoVFWMpm83MbdcvvtBhN/qumWlSMUahCxcnhPT7e8qoF2bc2ufeArgQAho4pVKdbWtr73qqBY9HDRoUZxnTaG/X3Df3sQpVVccDOv1hitQmiGsTZSJa8M5tC+8dOcGxnQfgGAIaOGI2m4qLC82XNyfqbEtOEh4eERgYRAjZumnH+2/+n92VhJTtj1HkwmhZU6Q4Izruz/uf8hd3+sYjAA+CgAYu1NXVXrrUcYDCzt7IgYFB4eERhJANH335yftfmZluHh5afoyhqNohfup4SX+F8viyVZH+gQ7uOoDrIKDB6VSqc21trI3orrqDpmlq8OBEmqY/3/jDB69v6OFr2Jr6iarS5JaXlMxMTN12/1OO6jCAm0BAg3OVlBQbDF1t7xke3i8wMLCsSHXfHY8ZdD3a1ag9SHBugv/lH1f4/9+dSycPSnRAXwHcDAIanMVsNhUVFTGdj1SIxeLY2MEmk3nujMUXzl/s0XfSpDJd0RomJITcOWL8GzcvUPphw2XwWghocAqDwVBaWsQwnU64iIyM9PdXlhSV3zv30b+eHHY1t8MkpMon+utlfD+B4IVJNz028Xqp0M4L+gC8CQIaHE+n05SVlRFC7OYtRfHi4+Npmv/D1l2vvPjBVU12jibNEeILaTKGIiEyxdb7Hk+PGuyUTgO4HwQ0OFhzc1Nl5fnOWqVSycCBgwgh77+56fNPtna1hwZFGqLE1SNkfIr3yZ0P3p6azud191pVAO+CgAZHunjxQn19pwu4AwKU/fpFMgxzx6yHVecqu0jn+mjJxeFSQkhCSMSJJ1/BnnPgmxDQ4DCVlZXNzfa3DKUoEhERGRCgLCtVLZi9zGzzLm0rnZx/dqo/Q1E8ivpiwaOzU0Y5rb8A7g4BDY6hVtd2ls6EULGxg0Ui0TP/8/L+PYc7/QqKnJsY0K7kE0JCZPKCZ9+WCARO6SuAh0BAgwNotZqamk53MoqPjzcazFPS7mhrabu6pcOsDZo6c2Mgw6MoQtbdseSekZlO7C6Ah0BAw7UymYzl5WftNlEUiY8fsuvH/f9a8aa95YGX09kooopnBRFCFGJJ0XPvyEWYPwdACAIarl1JSfFfi7OvmsdM0/zY2MG3TH2g5pK6iynOrWGCinH+hDArps1ZiX1BATpAQEPf/bUaxVpxJYIFAqFCqpw4fK7JZGY1daQaH9AWwufTdNGKt8Pk/s7tLoCn4bm6A+CpDAZ9SUmR3X2NeDzqj2NFN0+9/690tvfjEt6Z2UFtIXxCSO7TryGdAWzhDhr6qLS0pLOmV174+Ey+/VFpi4tJfvXxfoQQiiJbFjwWpQx2fP8APB8CGvri7NnSzjYFff+1zV2lM0XKpih1CpoQEiELyHv2DTHm0gF0AgENvdbc3KjVau21MD/9sP/k0fzOfrCln+h8upwQIhWJ9i97MTE00ml9BPAGCGjoHYYxX7hwwW7Tm2s25P1pf9yDoci5zABNEJ8Qsih98vu3LnJiFwG8BQIaeqesrMzu4MbK5W9Xquzv6WymqaIbgxgeIYS8N3fRA2MmO7ODAN4DAQ29oFbX6vU62/qd3+3tLJ01Sn75pABCSKBEunPJs8P6DXRuFwG8CAIaeqq9vdXueu56ddO2L3fb/ZGKcf6tYQJCyGs3L1g6/jrn9g/A6yCgoadUqgrbSoZhVv/9fbvHX0yRtYYJKEK+uOfRW1JGO7l3AF4IAQ09UlJi/+2CWzZsb2pssa2vGi5rjBbLRZLSF97Bu6kA+gYBDd0rLy8zGAy29Zeq6/b8fIRVyRCiygxoD+Y/OemmNbPmcdJBAO+EgIZuqNVqjUZjW28wGJ9/7A1WpV5Kl85QBvnJzzz9SpCfnJMOAngtBDR0xWQy1dTYn56x8vG3WVttXEr2q4/z+/XhleOi4zjpHYCXQ0BDV6qqKu3W7/x2z6WLakKIdR/RqmFSkhRUsnxtuELJYQcBvBkCGjql02lbWuw8ADyy/9S2r3796xNFCKkaLp92y6RNdz3CYe8AvB8CGjp19mypbWVdbeP6d7/uWNOS6F/6+SbseQTgcAhosK+yUmV3u7qVT7zd8SM/xO/Mf77kqE8APgYb9oMdWq2mudnO4Mah309pNFf2saNo3qGsbzjsF4BvQUADm9FoLC8vs6036A2fvHfV4MZnP7xDUfbfZQUA1w4BDWwlJcV2BzdYS7qHLJk8JH4QR30C8EmuCei6ujqXnBe6VVlZYXdJ97EDOZUVV3ZKahmi+N/Hl3PYLwBfxEVA19TULF68OD09fcWKFWq1euTIkeHh4QkJCYWFhRycHXqOYZjm5mbber3esP69ry4fQ4g2SHDym00iPp4wAzgXFwG9ePHi5ubmZ599tqCgICkp6cEHH9TpdI888shjjz3Gwdmh5yor7SxLYRjmmUf+bTZdHvXQBgp2/bzRTyDktmsAvoiLm6CsrKyKigp/f//09PT4+PgHHniAx+MtWbJkzZo1HJwdeshkMrW0NNnW/2vFB00Nl2d0NAyS7Prsw0j/QG67BuCjuLiDDg0NPXDgACFk//79Wq327NmzhJDi4uLAQPx37kZUqnO2laeOF5wtuXxbXZsiP/7dZ4NDwjjtFoAP4+IO+vXXX7/jjjv8/f0Zhvnwww9vu+22GTNm/Pjjj08++SQHZ4eeaG9v1WrtbFn3yftbLYWmKMn2919XSmTc9gvAp3ER0HPmzDl//nxpaWlKSopUKk1NTd2zZ8+HH344c+ZMDs4OPVFZed62MudEQXublhCik9HBNw1JjRjAeb8AfBpHD+KDgoKCgoIs5czMzMzMzMbGxt27d3eR0e3t7Rcvsje6bGpqsvtKabgWTU0NRqOJVWkymT94fQshhOGR8qkBpQ887YquAfg0l82UKigomDNnjt2d4C327t27efNmVuWff/45cCBeC+1Y5qqqC7a1P/2w12AwEkLKpgUlR0aFyvw57xiAr3NZQGdkZHSRzoSQm2666aabbmJVLl++vLq62pn98jllZaW2f5No2jXff/krQ0hjrFgvpX5c/Kwrugbg67heSdjS0oIxCvfR3Nyk0+lt6z9+52tCiJlPVafKHkifHCTDy6sAXICLgNZoNGvXrk1ISBCLxQqFQigUxsfHr1mzRq+3Ew3AIcbu4EbtpcacE4WEIqXXB0mFonfm3s95xwCAEG4CetmyZdnZ2Rs3bqyurjYYDDU1NZs2bSooKFi6dCkHZ4fO1NbWms3sbTfMZrLm2fcIIWVTlSY+2f8/q3nYrw7ARbgYg965c2dFRYVEIrF8VCqV48aNS09Pj42N5eDs0JnaWrVt5dZNP7Y0t9XHinVy+qPb/5YYGsF9xwDAgos76AEDBuzevZtVmZWVFRwczMHZwa7q6kpC2LfPzU0tu3ce0PnzL6bKJg5KvG/0BJf0DQAsuLiDXrdu3bx581atWpWUlCSXy1tbWwsLC9Vq9fbt2zk4O9gyGo0NDY229auf/cgkoMqmBEwZnLzzb3/nvmMA0BEXAT1mzJjy8vJ9+/apVKq6ujqlUrlkyZJJkybRNM3B2cFWaWmJbWVVZa1aXV98c/B9aRM/um0x970CABaO5kHTND19+nRuzgVd02o1ZjN73SAh5IPXNp/L8J+bmo50BnATeOWVz8nPL7CtLMgrPalsf/i2OZsXLOO+SwBgFwLat/x2+IhEImJVmoymf/3vlkWLb/vXrDtc0isAsAsB7UNMJlOoXGpbv2Pbb9Kbhqy5HukM4F4Q0D7k88+30TT7n7hWq3v/4qndS59zSZcAoAsIaF+x778H09OTbOtf2/jN9ytWykUS7rsEAF1DQPsEs5nJy7PzbLDmUv0/Vz2WEZPAfZcAoFsIaJ/wj7fXTZ6RzqpkGMYskSSGRrqkSwDQLQS09/vXf7cFGBg+nz3nvaBcNXVsmku6BAA9gYD2cofLiz//cOusm+zsqhE3PJn7/gBAzyGgvdxdb7zy+v0LeTz2P+izavXw/jEu6RIA9BAC2pttPP77kNOGuAQ7b3Gcms4ekgYAd4OA9mb/fuuTZcsX2NabTIyf2I/7/gBAryCgvdbP+X8EnWkdFD/Atikpyc6EaABwNwhor/X3J1+Zf++NtvVKZRA2egXwCAho7/RDznFZpXbydWNY9RRFwsPDXdIlAOgtBLR3WrnijVmzJ/P57Dtlf38lhZfAArgrdW0AACAASURBVHgIBLQXemrH536q9ptvn2zbhNtnAA+CgPY2l1qatq//z/Trx4lt9n0WiyU8HkafATwGAtrbTP/oX6GVprvuv8m2KSoqivv+AECfIaC9yseHfm06UTF1xhjb0WehUMjnC1zSKwDoGwS092jWap7e+UVofvsNcybZtkZG9ue+SwBwLRDQ3uPwuaKoIy1hwcqAQAWrSSgUSiRYOgjgYRDQ3uOj7Ttkl3RLn77btmngQOyLBOB5ENBegmGY819lE4qKGcxe283n8wUCjD4DeB4EtJd4fN1HAq05Zehg26aQkBDu+wMA1w4B7Q10RsPhdbsIIQuWzGY18XhEqQxyRacA4FohoL3BY8+8zDMwgUH+EZGhrCaZzN8lXQKAa4eA9ngGg/HUrmxCyD1L5ti2Bgbi9hnAUyGgPd6c+csshdSR8awmiqL9/DC7DsBTIaA9W11dQ21BJSEkY/JI2/d2BwcHu6JTAOAYCGjPdt/8JwlDCCEZE0eymiiKQkADeDQEtAervnCp9nwNIYRHUwnJ0axWqVSKrZ8BPBoC2oP9e/U6S2Hy9DG2S1FCQsI47xEAOBIC2oMdzsomhBDC3L7gettWiUTCcX8AwLEQ0J5qz+6DlsKQoYP8ZOwslkplnPcIABwMAe2p3n1jo6Ww8MFbbVsjIvpx2x0AcDwEtEeqU9dfOn+JECIUC8Ij2VM1BAK+UCh0Rb8AwJEQ0B5p1dNvWArXz55ICHuqBm6fAbwDAtrz1Ksbs4/mWMqzbpnIaqVpWiZjb9gPAJ4IAe15Vq94y1JISI6VSMSsVoUigPMeAYBTIKA9T84fpy2FuxfdbNsaGordnwG8BALaw+T8UdDeriOE8GhqYCx7rFkoFNA0e0cOAPBQCGgP896bn1oKM29ijz4TQsLDI7ntDgA4EQLak2g1uryTBZbyDXPYAc3j8WQyrE8B8B4IaE+yZMEzlkJUdITcn53FCgVengLgVRDQHqOy4mJRQZmlvOBv7HcPEkJCQ7E7EoBXcUFANzQ0tLS0cH9eT/fea/9nLSckxbJaBQKh7Yb9AODRuAjokpKSKVOmVFRUqFSqjIyMkJCQwMDAqVOnVlZWcnB272A2M8cO/2Epz50/3faAsDDcPgN4Gy4CeuHChRMnTuzXr99jjz02derUtra21tbWzMzMhx9+mIOze4fiM2Uay+w6Hu/6WybZtPMwAA3gfbj4o7iwsPD333/n8/m5ublbtmwRiUSEkH/+85+46eu5Tz780lKIHxItErE3QpLLMXkDwAtxcQc9duzYbdu2EUIyMjL27Nljqfztt9/69+/Pwdm9w4G9xyyFW+++zrYVuyMBeCUu7qA/+eSTWbNmvfnmm+Hh4bfffvv06dPNZnNBQcH333/Pwdm9wOGsE5YCj6bih8SwWoVCMR4PAnglLv7DjoyMzMnJ2b9/f0lJycSJE/39/aOjo6dNm4ZY6aG3/v2JpXDjrVNtWyMiwrntDgBwhKOIpChq8uTJkydPttY0NjYeO3Zs5syZ3HTAc1VfuFRRdnm6y9SZY1mtPB7l5yflvFMAwAWX3cMWFBTMmTNHo9F0dsC2bds+/vhjVmVxcfHgwYOd3DX38vT/vGwpSGV+ykD2Rs9SqZyi2Bv2A4B3cFlAZ2RkdJHOhJDbbrvttttuY1UuX768urramf1yL1qNtuTM5dWDM2/OtD0gPBzjGwBei+uVhC0tLQzDcHxSz/Xs469ay5NnpLNaBQKBQIB3DwJ4LS4CWqPRrF27NiEhQSwWKxQKoVAYHx+/Zs0avV7Pwdk9l1aryz5yylIOjQhWBMhZBwQEKDnvFABwh4uAXrZsWXZ29saNG6urqw0GQ01NzaZNmwoKCpYuXcrB2T3Xpx9vNRqMlvL8+24ghP2XR2BgEOedAgDucDEGvXPnzoqKColEYvmoVCrHjRuXnp4eG8ve8Qc6+mbLTms5PimG9fZugUBI0zTnnQIA7nBxBz1gwIDdu3ezKrOysoKDgzk4u4c6nVvU1tJuKQco5TK5H+sA/PYAvB4Xd9Dr1q2bN2/eqlWrkpKS5HJ5a2trYWGhWq3evn07B2f3UBs/3Gotj8kcbnuAUhnIYXcAwAW4COgxY8aUl5fv27dPpVLV1dUplcolS5ZMmjQJf6F34WT2n9bydTdPYLX6+bFvqAHA+3A0D5qm6enT7exiDHYdOXDSsrkoISQoWBkUzN5KVCZjr1gBAO+DV165o//9YIu1fN9Dc2wPUCoDOOwOALgGAtrtMAxTVqyylPl8etioRNYBfD6fprHPFID3Q0C7nRPHcrQaraU8e56d7euCg0O57REAuAYC2u189NZma3nSjDGsVoriBQZi/gaAT0BAuxej0ZifW2QpBygV/jbLu6VSbC4K4CsQ0O4la99xazktI8X2ALwnDMB3IKDdy1efXn4NGENIWkYqq1UoFPJ4mDwO4CsQ0G7EbDIXF56zlHmUndcP+vtjdh2AD0FAu5GiwjJN++WXGIzJHGZ7QEAAAhrAh3QV0AaDgbN+ACFk21c/W8vTZ2WwWmmaxvb8AD6Fvd5Br9dv3Lhxz549ubm5NTU1Eolk6NChaWlpjz76KN6u5GwHf8+2lkPD2Xs9S6UybrsDAC521R305s2bZ86c2dra+tRTTx05cqSxsbGsrOyNN95ISkpatGjRc88919zc7KqOer262vp6dYOlHBgcoAhgx3FQEPYXBfAtV91BS6XSPXv2dNxkzs/PLzU1NTU1dcGCBX/88UdNTY1CgW16nGLXjt+t5eGjh7BaaZq2vvEAAHzEVXfQt912G03TP/74I+ugzz//nBAycuTIwYMHc9c1H7P18yu7Y0+ekcZqFYuxvyiAz7HzkPDFF1+8/fbbL1y4QAjJzc2dMGHCBx98wHnHfEtjQ/Oli2pLmc+no2IiWQfIZAhoAJ9jJ6Czs7OnTp06fvz4+fPn33DDDYsXLz569Cj3PfMp33xx5fWDo8cOtT1ALmdvCQ0AXs9OQPN4PKVSyTCMWq0ODAwcPHgwRVG2h4ED/bx9r7U80WZ8g6ZpoRAT7AB8jp1thSdPnmwymX766aeUlJSsrKyHHnooLS3ts88+47xvvsJgMFZVXrR+TBrKHujHg1kA32TnDnrhwoVZWVkpKSmEkIkTJ546dQrPBp2q5Ey5tTxgoJ3J5ng/LIBvuiqg9+7dSwhZtGhRxzENoVC4cuVKQkhpaalKpeK4f77gg7c+tZbtvsBbLMYEOwBfdFVA5+XlXXfddZ9++mlRUZFlnTfDMJWVlbt27Zo3b95bb70lk2Exm+MV/LUBNCEkfTx7Cw65nL0lNAD4iKvGoB9//PHFixe/9957jzzySEFBgUgkamtri4mJSUtLe/XVVwcNGuSqXnqx0qLy9vbLL7iS+IlDw9mjGQoFNkgC8FHsh4Qymez5559//vnnCSF1dXVKpZLHw453TvTLj/ut5SEpdv4XiL9aAHxWVy+HDgpi79cDDrf3vwet5ckz0lmtQqGw48p7APApdgJ6z549q1evbmxs7FiZl5fHVZd8iF6nr1RVWz+mjkokhBDCEHL5IS0GoAF8mZ2AXrx48ZIlS2699VbcuznbeVW1NY7jrrw/5coUGiwgBPBldgJao9G88MILWD3Ige+3/myN4wlTR7FaKYrn54ctOAB8l50HgNOnT1+/fj3DMNz3xqcwDLNz23+tH1OGxbEOkMvxeBDAp9kJ6IqKikceeSQ4ODgxMXHIX7jvmdc7cTRHq9VbylKZX2Awezqdv7+S804BgBuxM8Sxfv167vvhgzZ+vNVanjCVvUESwQQ7AJ93VUAHBwer1erk5GRCyJdffnn33Xe7qFfez2wynzldbP04amwy6wChUITHAAA+7qohjqamJmt54cKFnHfGh5SXVbS3aa0f4xIHsg7ADnYAgFWCrvH7b1fegRAWYX0b7JUHs4GB2MEOwNchoF3jxNEca3nsROsOdpfHNPh8Pp8v4LxTAOBerhqDZhgmJyfHtkwIGTaMvcsaXIuTx3Kt5YyJI1it2F8UAAgroJVK5fTp023LhJDa2lpO++XVTp3It5Z5PKrDEMdlmL8BAIQV0Ehhbuz79dBfRSZm8MCOm29YKJWYAQ0AGIN2hexDf/5VpKbfkMFKZ4FAQFH45wIACGjO6fWGS9U11o/Jw9jve5RIsP8GABCCgOZeYX5pS5vGUpYrpAp/9nAzxjcAwAIBzbUff/jVOqIxcbqdFd5+flIu+wMAbgsBzbXDWSet5eRU9vgGVngDgBUCmlM6rb7motr6MSqmP+uAwECMbwDAZQhoTuXlFFrLAwaGy+TsBSl4hzcAWCGgObX3F+sMaHL9LZNsD+Dzu3qNLwD4FAQ0p04ev7J6PmEIewc7qRSPBwHgCtcE9JEjR1xyXtcym5nysvOWslAoCA4LYh2gVGIHOwC4wjUBfd1117nkvK5VXHjWWo6KjmC10jQll2MPaAC4gouAfuKJJyRXa21ttRQ4OLv7OLj3mLWcmMKeYMfnY4IdAFyFi0dSa9eu1Wg0+fn569atCw0NJYTExsaePXu22x/0MgWnS6zltAz2O64wAA0ALFzcQfv5+a1fv/7vf//73XfffeLEibCwMB6PFxYWFhYWxsHZ3UdhQam1PDCWPQMaAQ0ALNxN6po9e3Z6evqiRYt27tzJ2UndR2tLm7qm3lKWK+xksVSKPaAB4CqcPiQMDw//+eefhw4dOmfOHC7P6w4O7jtuLQ8ZOojVKhAIeTxMeQSAq3AdChRFLVu27PPPP29sbNy9ezfHZ3eh//6UZS2PHjuU1SoWi7ntDgB4AJetWysoKJgzZ45Go+nsgNzc3KNHj7Iq8/LyPHTuR1npOWs5KrofqxV7QAOALZcFdEZGRhfpTAgRi8W2OyOLxWJPHArQanXVFy5v0i+WiMIj8RJCAOge1wHd0tIik8l6MuE3Pj4+Pj6eVXn48OHq6mrndM2Jck5eeUvs4AT2Cm+a5olEIm57BAAegIu7UY1Gs3bt2oSEBLFYrFAohEJhfHz8mjVr9Ho9B2d3B7//emVp+/DRiaxWgUCIJSoAYIuLgF62bFl2dvbGjRurq6sNBkNNTc2mTZsKCgqWLl3KwdndQceATh7G/rMAr1ABALu4GOLYuXNnRUWF9eGeUqkcN25cenp6bGwsB2d3ucaG5rq6BkuZR1MRkSGsA7AFBwDYxcUd9IABA2xn1GVlZQUHs5+VeaU/juUSwljKdpcLYg0hANjFxR30unXr5s2bt2rVqqSkJLlc3traWlhYqFart2/fzsHZXa7gdAkhl4eYE1NiWK0iEWZAA4B9XAT0mDFjysvL9+3bp1Kp6urqlErlkiVLJk2aRNM0B2d3udxTZ6zluMRoViuWqABAZziaZkfT9PTp07k5l7v5s8s5dh667gYAOOB5iz48S1mxylpmCImNG8A6AFM4AKAzCGjnOv3nldd4h4bZeaOVUCjksDsA4EkQ0M516sRpa7n/QPZrrgQCgSeuXAcAbiAdnKvyQpW1HJ8UzWrFHkkA0AUEtBOZzYyq9IL1Y9rYFNYBfn4IaADoFALaiSrOVTY1NVvKQhE/OJQ9Bo1J0ADQBQS0E/2264B1iYq/v5313LiDBoAuIKCd6PiRXGvZdgBaLBZjEzsA6AIC2olOncizlqNiIlmteEIIAF1DQDtL5fmrXiwQHWv7miusIQSAriCgneWP43kdPyYks/dWxWuuAKBrCGhnObQ/21r2k9q5WcYSFQDoGjLCWWouqq3lhCT2LqMCAZ/H84nN/ACgzxDQzlJX22Atj0hPYrWKRBiABoBuIKCdorqq5mJ1rfVj7GD2JnZiMQIaALqBgHaKvbsPWct8Pt1/YDjrADwhBIBuIaCdIq/DW1QSkuy8Gxdz7ACgWwhopygsKLWWE1IsAc1Ya4RCIdYQAkC3ENCOp9fpqyovWT9GRVu2gb6SyHiLCgD0BALa8Vpa2jt+DO8XzDogIEDJYXcAwFMhoB0v91SBtSwQ8MP7hbAOwJu8AaAnENCOt+O7/1rLQ4YOsj0AawgBoCeQFI6nKq+0loenDWG1Yg9oAOghBLSD6fWGyoor+9jFJ7IXeeMJIQD0EALawfJziqxlgYCPJSoA0GcIaAc7W6qylgfGsPeAJljkDQA9hoB2sD9OnLaWh6fZ7pEkwhNCAOghhIWD/Xki37pkMDY+itWKFd4A0HMIaEdqbW6rvVRnXTKYZDPHTqHw57hLAOC5ENCO9J9vd1vLYrHI9gBM4QCAnkNAO9K5sxXWctyQaFarUCjEADQA9BzywpHy80qs5bGZw1itAoGA2+4AgGdDQDtSWfE5a3mAzRw7jG8AQK8goB3mwvmLHT/+tcvoFVIpAhoAegEB7TDHD/9pLQcGB9geIJFgFw4A6AUEtMMc2n/cWh5hs0RFLBbhLSoA0CsIaIepulBjLcclspeo0DSeEAJA7yCgHcNgMJYWlVs/BocGsg6QSrFHEgD0DgLaMc6Vne/4cXDCQNYBCoWCw+4AgDdAQDvG6ZxCazk0LMj2AKFQyGF3AMAbIKAdo7jwyvhG0jD2FhxCoZ1l3wAAXUNAO0ZlRZW1PDZzBCGEEOuudkQkwu0zAPQaAtoxjh06ZS2HhFmeEF6ZVCeTyTnvEQB4PAS0A7S3a6xlhb8sKIS9SkUuxxQOAOg1BLQDnDyWay0Pstmkn6ZpPh9DHADQawhoB/jzRL61LLQZbvbzwwpvAOgLBLQDHDmWYy0np7KncMhkmAENAH3BXUA3Njayapqamjg7u1OVni61luOTYlmtfD7NbXcAwEtwEdB5eXnJyclBQUGDBg3avn27pVKn0wUE2NnyzeM0NbZYywIBP7xfMOsAkQiToAGgL7gI6IcffnjRokVarXbDhg0PP/zwoUOHODgpZ0o7bNJvM3+DIVilAgB9xefgHPn5+QcOHODxeJMnT16/fv2SJUtycnK6/zEPUVxy7q8i4x/Amu9M4QkhAPQZF3fQAwcO/O233yzl2bNnp6amLl++nIPzcuPIcesSFSpj0siOCwgJIRKJhPsuAYB34CKgX3311blz506YMKG+vp4Q8vHHH2dnZ48fP56DU3PgouqStZyUOrjjAkJCiECAGdAA0EdcDHHMmjWrsLDw4MGDlsdlAQEBBw8e3L59+4kTJzg4u1MxDHOu8Jz141+LvK/AHTQA9BkXAU0I6d+///z5860fBQLB9OnT5fKudqgwGo0tLS2sSp1O55T+9VVV5ZW3qMgVdt4Ji/cQAkCfcRTQtgoKCubMmaPRaDo7YMeOHevWrWNVFhcXx8XFOblrvfBH9pVF3gk2M6DFYtw+A0DfuSygMzIyukhnQsitt9566623siqXL19eXV3tzH71TmHRWUuBISTMZga0WIwJdgDQd1wv9W5paWEYpvvjPMTxE3mWaRuUnddcMRjfAIBrwUVAazSatWvXJiQkiMVihUIhFArj4+PXrFmj1+s5OLtTqfLLrdM2/KSs+2UKS1QA4FpwEdDLli3Lzs7euHFjdXW1wWCoqanZtGlTQUHB0qVLOTi7UzF+V8aIAgLYmyJhkTcAXAsuxqB37txZUVFhnXCmVCrHjRuXnp4eG8t+quZZSorLSbvRUlb4y1hj0CKRiM932RA/AHgBLu6gBwwYsHv3blZlVlZWcDD7qZpn2XbokHVRSmRUGKsV6QwA14iLEFm3bt28efNWrVqVlJQkl8tbW1sLCwvVarV1ZzsPtSc721pOSIphtXrRo1AAcA0uAnrMmDHl5eX79u1TqVR1dXVKpXLJkiWTJk2iac/eKLntfL21nJDMHq6RyfAeQgC4Jhz9GU7T9PTp07k5F2f0TVrr/2GGpLBfpIJF3gBwjfDKqz6qamqga9otZYHAzv/n/PzsrPwGAOg5BHQf/dJhpyfbJ4RCoYiiKAIAcA0Q0H10plJlLadlDGO1CoUCbrsDAF4IAd1HecfzreVhI+NZrSKRmNvuAIAXQkD3hclsvlBcaSnTNN1/YATrAOzCAQDXDgHdF8dUJXTD5Z2pBwxkD0ATQmQyPCEEgGuFgO6LE2dLha0mS9l2l1GhUMTjefYUbwBwBwjovvj94DFrOS4xmtWKN3kDgEMgoPtCXddoLYeGs++gsUQFABwCAd0XzWVXXkUYFh7EasUdNAA4BAK618431pnqLr+sSyqT2I5BY44dADgEArrX1M1NksbL20AnJrO34AAAcBQEdK/tPHLUWh42OoHVijd5A4CjIKB7rfRshbVsO4VDoWC/+AoAoG8Q0L12If/yLhx8Ph0RGcpqxSZ2AOAoCOheUxddsBQGxkbatmKOHQA4CgK6dwxGk6jl8hrC1JHsAWiRCLuMAoDDIKB753DZGWs5KqYfqxV7JAGAAyGge2f77v3WclQMe4hDJpNz2x0A8GYI6N75M7/I8rZusVgUFOzPapVK8YQQABwGAd072op6yxhzWAR7hTchxNPfUw4AbgUB3QtmhtE2tlnKQ0ewnxBi/gYAOBYCuhfOXDjPb9BbymMnDGe1YgY0ADgWAroXiorOWsv9B4Zf3chgCgcAOBYCuhdycgotBWWQ7XpuSiaTcdwfAPBuCOheOHT4D0thZHoyq0kg4PN4+GUCgCMhU3qhSXV5n/5JM8awmsRijG8AgIMhoHuqsbWVajEwhPD5dFR0BKsVA9AA4HAI6J5SNzYQQihCBiUMsG3198cuowDgYAjontr43U5LYdTYVFYTTdMCgZDzHgGAl0NA99TBo389IUwbwmqSy7EFBwA4HgK6p/RVLYQQHk0FhwaymhQK9qYcAADXDgHdU+aaVkJIVLSdTfoxAxoAnAEB3SO5JaWWwqTpaawmoVBECDbpBwDHQ0D3yP7jpyyFYaMTWU3YYhQAnAQB3SP1ZRcJIVKZX2BQAKsJm/QDgJMgoHuk4twFQkjK8DjbJtxBA4CTIKB7pLTkHCFkhM0WHEKhCFtwAICTIFy6p9cbNLUthJDk1EGsJoUCCwgBwFkQ0N2rb2gihASFKuUK9nQ6BDQAOA8CununcgoIIROmjrZtEovxmisAcBYEdPf+OJpLCJkwZRSrHi8hBACnQkB379ChExKpOChEyar398cKbwBwIgR099Sq2qHD4m3rAwLYm3IAADgQArobLc2thJAxmcNY9UKhEBPsAMCpEDHdyPuzkBCSkBzDqscWowDgbFwHdEtLC8MwHJ/0WhRfrAqPCJLJ2csF/f3Za74BAByLi4DWaDRr165NSEgQi8UKhUIoFMbHx69Zs0av13Nw9mt07Pifc+6cYVuPCXYA4GxcBPSyZcuys7M3btxYXV1tMBhqamo2bdpUUFCwdOlSDs5+jYryy5KGsbfgwAQ7AOAAn4Nz7Ny5s6KiwhpqSqVy3Lhx6enpsbGxHJz9GtFak8KfvYAQ8zcAgANc3EEPGDBg9+7drMqsrKzg4GAOzn4t6tQNN8ya0KHi8ui5UsmeEw0A4HBc3EGvW7du3rx5q1atSkpKksvlra2thYWFarV6+/btHJz9WtQ2N02a1nGFN0UIEYvFruoPAPgULgJ6zJgx5eXl+/btU6lUdXV1SqVyyZIlkyZNommag7Nfi33HT05KTWBVyuVYQAgAXOAioAkhNE1Pnz69Y01jY+OxY8dmzpzJTQf6hm802lYGBmIAGgC4QLlqVvLhw4enTZum0Wg6O2DHjh1ffPEFqzInJ+fFF9eIREKpVKI1mNs0Oj+RQCK6/L8Zh9fI/ER8Pk0ICQjwF4lEer1Bra6zzA6MjY1NTEwsLy8/c+aM5adQgxrU+GDNiRMnRo4cOWgQe7N4h3BZQHfLYDC0trayKr///vv65tZb5t7CMIxIJKJ5PJPZrNPpLK0Oryk5Wx4RICeEUBSPogjDMEYeHRPZnxAilUppmjaZTG1tbZafQg1qUOODNU899dSjjz7qJQHd0tIik8koiurbj3/zzTdqtZrLCdR7jx0Nl1+eZlev12cOH8nZqQHA/f3tb3977rnnnBTQXIxBazSat99+e9OmTSqVSqfT8fn8mJiYe+65Z8WKFUKhkIMOXIupY8aazWaDwSAQCLA7EgBwCSsJu8fj8UQivBwWALiGlYQAAG4KKwkBANwUVhICALgprCQEAHBTLltJCAAAXcPMBAAAN4WABgBwUwhoAAA3xdEYtKP4+/uvXbv2hx9+4PKk1dXVFRUVAoGAy5O6nMlkMhqNIpHI1R3hlNls1uv1vrblN8MwOp3O167aYuTIa9284fz58857B577bpbkPrZt21ZVVfXoo4+6uiOc+v3333///fd//vOfru4Ip3Jzczds2PDuu++6uiOcunDhwrPPPmu7eaR302q1c+fO3bVrl6s70hUMcQAAuCkENACAm0JAAwC4KQQ0AICbQkADALgpBHT3aJr2wW1DfPOqeTyeD2787ZtXTVGU+/8bjml23TMYDCaTyddmiZrNZo1GI5VKXd0RTjEM09bWJpPJXN0RrrW0tMjlclf3gmvuf9UIaAAAN+Vzf9cAAHgKBDQAgJtCQAMAuCkENACAm0JAAwC4KQQ0AICbQkADALgpBHQ3jh8/PmrUqJCQkMWLF2s0Gld3x/F27do1fPhwpVKZmZl5+vRpS6Xdq/a+X0V7e/sTTzxh/ej1V11VVXXDDTcEBQVlZmaWlJRYKr3+qrdv356QkBAQEHDLLbfU1tZaKj3mqhnonNFoHDBgwPfff9/e3n7LLbe8/PLLru6Rg50/f14qle7evdtoNL722mtJSUlMJ1ftlb+KBx98UCwWW8pef9Vms3no0KFbtmzR6XT/+Mc/Jk+ezPjAVavVaolEsmPHjoaGhttvv/2BBx5gPOqqEdBd+fXXX1NSUizlgwcPxsfHu7Y/DvfVV19lZmZaynq9nsfjNTQ02L1q7/tVbN++ffz48daA9vqr3rdv36hRoyxlg8GgUqkYH7jqo0ePhoeHW8rfffddWloa41FXjSGOKxAUNAAABAxJREFUrqhUquTkZEs5KSlJpVK5tj8Od/PNN2/fvt1SPnToUHR0dEBAgN2r9rJfxcWLF1944YV169ZZa7z+qktKSiIiIhYsWBAbG3v77bdbKr3+qpOTkxmG2bJlS2Vl5aZNm6ZOnUo86qoR0F2pq6uz7qUil8t1Ol1LS4tru+RYUqk0MDCQELJjx4677rrrrbfeIp1ctTf9KhiGWbRo0auvvhoaGmqt9Pqrbmho2LVr1+23356Tk5OcnDx//nziA1ctk8lWrVp1zz33xMXFHT169JlnniEeddUI6K4olcrW1lZLuaWlhc/ne98+Z01NTXfdddezzz773Xff3XLLLaSTq/amX8VHH300cODAG2+8sWOl11+1v7//tGnT5s6dK5fLV61adfz48YaGBq+/6oMHD37wwQcFBQVNTU2rV6+eMWMG8ah/1gjorsTGxhYVFVnKhYWFMTExFEW5tkuOZTAYZs6cGRgYmJubO378eEul3av2pl/F0aNHt2zZIpfLY2NjtVqtXC4/cuSI1191dHR0x48URfH5fK+/6j179lx33XVDhgwRCoUPP/xwXl5eXV2dJ121qwfB3ZrRaIyMjNy3b5/RaJw/f/7q1atd3SMH+/bbb0eNGqXtgOnkqr3yV3Hx4sWOszi8+6r1en1kZOSePXtMJtOqVaumTJnC+MBV79q1a9CgQadOnWpra3v77bctj/486KoR0N04ceLEiBEjoqKiFi1apNPpXN0dB1uxYgXrf9iNjY1MJ1ftfb+KjgHN+MBVHz16dPjw4cHBwbNmzTp//ryl0uuv+oMPPoiJiQkICJg2bdqZM2cslZ5y1diwHwDATWEMGgDATSGgAQDcFAIaAMBNIaABANwUAhoAwE0hoAEA3BQCGgDATSGgAQDcFAIaAMBNIaABANwUAhoAwE0hoAEA3BQCGgDATSGgAQDcFAIaAMBNIaABANwUAhoAwE0hoAEA3BQCGgDATSGgAQDcFAIaAMBNIaDBt3z66adTpkwhhJjN5lGjRv3000+u7hFApyiGYVzdBwDuMAwzceLEpUuXNjQ07N2797vvvnN1jwA6hYAGn3P69Okbb7yRoqhDhw5FRka6ujsAneK7ugMAXEtJSYmLiwsODkY6g5vDGDT4nKysrJqamoMHD+bm5rq6LwBdwRAH+Ba9Xj9ixIj169eXlpZu2LDhwIEDru4RQKdwBw2+5bXXXktLS8vMzFy4cKFer9+8ebOrewTQKdxBAwC4KdxBAwC4KQQ0AICbQkADALgpBDQAgJtCQAMAuCkENACAm0JAAwC4KQQ0AICbQkADALgpBDQAgJtCQAMAuCkENACAm0JAAwC4KQQ0AICbQkADALgpBDQAgJtCQAMAuCkENACAm/p/sLMeubQp+v4AAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "###----------------------------------------------\n",
    "## matching \n",
    "\n",
    "# filter Noncausal to match SNP centrality\n",
    "centr = NONLEAD_tabf0502$Dist_SNP_to_peak_center\n",
    "uloc = unique(LOCAL_tabf_2$Peak)\n",
    "ulocdis = LOCAL_tabf_2$Dist_SNP_to_peak_center[match(uloc,LOCAL_tabf_2$Peak)]\n",
    "print(length(uloc))\n",
    "min_v=c()\n",
    "for (i in 1:length(ulocdis)){\n",
    "  w=which.min(abs(centr-ulocdis))\n",
    "  min_v[i]=w\n",
    "  centr[w]=1000000 # set to large number to avoid resampling\n",
    "}\n",
    "\n",
    "NONLEAD_tabff = NONLEAD_tabf0502[min_v,]\n",
    "print(dim(NONLEAD_tabff))\n",
    "print(length(unique(NONLEAD_tabff$Peak)))\n",
    "\n",
    "plot.ecdf((CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center+1),col=enh_cols[3])\n",
    "plot.ecdf((LOCAL_tabf_2$Dist_SNP_to_peak_center+1),add=T,col=enh_cols[2])\n",
    "plot.ecdf((NONLEAD_tabff$Dist_SNP_to_peak_center+1),add=T,col=enh_cols[1])\n",
    "\n",
    "print(wilcox.test((LOCAL_tabf_2$Dist_SNP_to_peak_center),(NONLEAD_tabff$Dist_SNP_to_peak_center)))\n",
    "print(wilcox.test((LOCAL_tabf_2$Dist_SNP_to_peak_center),(CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center)))\n",
    "print(wilcox.test((CMLEADS_SNPinpeak_2$Dist_SNP_to_peak_center),(NONLEAD_tabff$Dist_SNP_to_peak_center)))\n",
    "\n",
    "GR_nonleads = make_GR_from_df(NONLEAD_tabff,c(\"Chr\",\"Pos_Left\",\"Pos_Right\"))\n",
    "print(length(GR_nonleads))\n",
    "print(dim(NONLEAD_tabff))\n",
    "## save R objects \n",
    "#save(Gaf_LEADS_SNPinpeak_2,Gaf_LOCAL_tabf_2,Gaf_NONLEAD_tabff,file=\"home/.../indiviudal_peak_tables_with_SNP_enhancers.RData\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1850240",
   "metadata": {},
   "source": [
    "### Combine everything to a complete table and add neighboring peaks of CM-dependents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "87007c0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "## First look independent of SNPs in GM for a general characterization\n",
    "vcm_Lids = unique(CMLEADS_SNPinpeak_2$Peak)\n",
    "\n",
    "LEDE_tab_vcm_1pair_only = CMLEADS_SNPinpeak_2[match(vcm_Lids,CMLEADS_SNPinpeak_2$Peak),]\n",
    "\n",
    "joint_tab = rbind(LEDE_tab_vcm_1pair_only,LOCAL_tabf_2,NONLEAD_tabff)\n",
    "\n",
    "# get Peak IDs for dependent peaks\n",
    "CMlead_dep_list_a = tapply(dag[,2],dag[,1],c)\n",
    "\n",
    "#local\n",
    "LOCAL_nei_list_a = lapply(LOCAL_tabf_2$Peak,function(u){\n",
    "  c((u-1),(u+1))})\n",
    "names(LOCAL_nei_list_a)=LOCAL_tabf_2$Peak\n",
    "\n",
    "#nonlead\n",
    "NONLEAD_nei_list_a = lapply(NONLEAD_tabff$Peak,function(u){\n",
    "  c((u-1),(u+1))})\n",
    "names(NONLEAD_nei_list_a)=NONLEAD_tabff$Peak\n",
    "\n",
    "## join all neighbor lists\n",
    "DE_list_a = c(CMlead_dep_list_a,LOCAL_nei_list_a,NONLEAD_nei_list_a)\n",
    "\n",
    "## combine lead and dep info\n",
    "LEDE_taba = add_dep_info(joint_tab,DE_list_a,vars,peaks)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "d38c3cec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Peak  Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7        7 chr1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222    213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.1  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.2  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287    277 chr1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296    286 chr1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "           Beta       SE  Log10_BF Chr Pos_Left Pos_Right Peak      GC\n",
      "7      0.886884 0.188153  3.085214   1   565005    565787    7 43.0396\n",
      "222    0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.1  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.2  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "287   -0.654595 0.049849 19.692796   1  3819984   3820189  277 58.7379\n",
      "296   -0.666789 0.041505 25.445000   1  4327332   4328105  286 53.8760\n",
      "      Peak_Height  P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7        13322207 0.985192 0.6912524  0.624720158        782      565396\n",
      "222          7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.1        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.2        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287          1990 1.000000 0.9377345  0.020502304        205     3820086\n",
      "296          9809 1.000000 0.9997858  0.004738315        773     4327718\n",
      "      Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr dep_Pos\n",
      "7                         110      0.1406650        8    chr1  569204\n",
      "222                       123      0.2026359      209    chr1 2508343\n",
      "222.1                     123      0.2026359      212    chr1 2508343\n",
      "222.2                     123      0.2026359      214    chr1 2508343\n",
      "287                        81      0.3951220      279    chr1 3820005\n",
      "296                        83      0.1073739      284    chr1 4327635\n",
      "         dep_RsID dep_Ref dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead  dep_Beta\n",
      "7     rs112660509       T       C 0.40125               1  0.4776132  0.251458\n",
      "222   rs200553981       G       C 0.35000               2  0.6725081 -0.134441\n",
      "222.1 rs200553981       G       C 0.35000               2  0.9016126 -0.205952\n",
      "222.2 rs200553981       G       C 0.35000               2  0.9806876  0.043439\n",
      "287    rs10909828       G       A 0.34500               2  0.9403693 -0.284635\n",
      "296     rs2101576       C       T 0.55000               2  0.9996411 -0.342060\n",
      "        dep_SE dep_Log10_BF dep_Chr dep_Pos_Left dep_Pos_Right dep_Peak  dep_GC\n",
      "7     0.140264   -0.6716076       1       569112        569721        8 44.0984\n",
      "222   0.032999    1.3310342       1      2495782       2496858      209 64.3454\n",
      "222.1 0.046273    2.0625997       1      2499680       2500108      212 66.4336\n",
      "222.2 0.013210   -0.1630566       1      2508970       2510224      214 62.0717\n",
      "287   0.051033    4.0437607       1      3826501       3826784      279 49.6479\n",
      "296   0.049869    6.4986698       1      4233102       4233752      284 52.8418\n",
      "      dep_Peak_Height dep_P_caQTL      dep_PMR dep_P_eQTL_Coloc dep_Peak_width\n",
      "7            29920971    0.015921 6.650092e-07      0.015463689            609\n",
      "222               402    0.032294 0.000000e+00      0.014451237           1076\n",
      "222.1             660    0.091908 0.000000e+00      0.042296138            428\n",
      "222.2            6161    0.012672 0.000000e+00      0.001904334           1254\n",
      "287               556    0.997449 3.472557e-03      0.172399409            283\n",
      "296               875    0.970727 0.000000e+00      0.006162762            650\n",
      "      dep_Peak_center dep_Dist_SNP_to_peak_center dep_SNP_centrality\n",
      "7              569416                         212          0.3481117\n",
      "222           2496320                       12023         11.1737918\n",
      "222.1         2499894                        8449         19.7406542\n",
      "222.2         2509597                        1254          1.0000000\n",
      "287           3826643                        6638         23.4558304\n",
      "296           4233427                       94208        144.9353846\n",
      "    Peak  Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7      7 chr1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287  277 chr1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296  286 chr1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "404  390 chr1 6837155  rs74050914   T   C 0.1350           1 0.8056400\n",
      "574  549 chr1 8960153  rs12034719   G   A 0.2550           1 0.9481063\n",
      "         Beta       SE  Log10_BF Chr Pos_Left Pos_Right Peak      GC\n",
      "7    0.886884 0.188153  3.085214   1   565005    565787    7 43.0396\n",
      "222  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "287 -0.654595 0.049849 19.692796   1  3819984   3820189  277 58.7379\n",
      "296 -0.666789 0.041505 25.445000   1  4327332   4328105  286 53.8760\n",
      "404 -0.528177 0.062199  9.967820   1  6836848   6837397  390 49.4545\n",
      "574 -0.409025 0.035001 16.526920   1  8959750   8960472  549 52.2822\n",
      "    Peak_Height  P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7      13322207 0.985192 0.6912524  0.624720158        782      565396\n",
      "222        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287        1990 1.000000 0.9377345  0.020502304        205     3820086\n",
      "296        9809 1.000000 0.9997858  0.004738315        773     4327718\n",
      "404        1152 1.000000 0.9589177  0.057274767        549     6837122\n",
      "574        2565 1.000000 0.6308161  0.035874555        722     8960111\n",
      "    Dist_SNP_to_peak_center SNP_centrality\n",
      "7                       110     0.14066496\n",
      "222                     123     0.20263591\n",
      "287                      81     0.39512195\n",
      "296                      83     0.10737387\n",
      "404                      33     0.06010929\n",
      "574                      42     0.05817175\n",
      "[1] 36119    50\n",
      "[1] 17287\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "print(head(LEDE_taba))\n",
    "print(head(joint_tab))\n",
    "print(dim(LEDE_taba))\n",
    "print(length(unique(LEDE_taba$Peak)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75d0abc1",
   "metadata": {},
   "source": [
    "## Add  additional information\n",
    "\n",
    "#### distance between peak pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "f31f4555",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "       CM     LOCAL no-effect \n",
      "     2051      7625      7611 \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "LEDE_taba$LD_peak_dist = abs(LEDE_taba$Peak_center - LEDE_taba$dep_Peak_center)\n",
    "\n",
    "vcm_Lids = CMLEADS_SNPinpeak_2$Peak\n",
    "loc_Lids = LOCAL_tabf_2$Peak\n",
    "nl_Lids = NONLEAD_tabff$Peak\n",
    "\n",
    "LEDE_taba$enh_subcategory = \"notassigned\"\n",
    "LEDE_taba$enh_subcategory[LEDE_taba$Peak %in%CMLEADS_SNPinpeak_2$Peak ] = \"CM\"\n",
    "LEDE_taba$enh_subcategory[LEDE_taba$Peak %in%LOCAL_tabf_2$Peak ] = \"LOCAL\"\n",
    "LEDE_taba$enh_subcategory[LEDE_taba$Peak %in%NONLEAD_tabff$Peak ] = \"no-effect\"\n",
    "\n",
    "print(table(LEDE_taba$enh_subcategory[!duplicated(LEDE_taba$Peak)]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87fe193a",
   "metadata": {},
   "source": [
    "### add annotation data PeakAnno & TAD information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "b1704958",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: snapshotDate(): 2021-05-18\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> preparing features information...\t\t 2023-04-05 03:31:44 PM \n",
      ">> identifying nearest features...\t\t 2023-04-05 03:31:44 PM \n",
      ">> calculating distance from peak to TSS...\t 2023-04-05 03:31:44 PM \n",
      ">> assigning genomic annotation...\t\t 2023-04-05 03:31:44 PM \n",
      ">> adding gene annotation...\t\t\t 2023-04-05 03:31:46 PM \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: 'select()' returned 1:many mapping between keys and columns\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> assigning chromosome lengths\t\t\t 2023-04-05 03:31:46 PM \n",
      ">> done...\t\t\t\t\t 2023-04-05 03:31:46 PM \n",
      ">> preparing features information...\t\t 2023-04-05 03:31:46 PM \n",
      ">> identifying nearest features...\t\t 2023-04-05 03:31:46 PM \n",
      ">> calculating distance from peak to TSS...\t 2023-04-05 03:31:47 PM \n",
      ">> assigning genomic annotation...\t\t 2023-04-05 03:31:47 PM \n",
      ">> adding gene annotation...\t\t\t 2023-04-05 03:31:48 PM \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: 'select()' returned 1:many mapping between keys and columns\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> assigning chromosome lengths\t\t\t 2023-04-05 03:31:49 PM \n",
      ">> done...\t\t\t\t\t 2023-04-05 03:31:49 PM \n",
      "\n",
      "        3'         5'     Distal Downstream       Exon     Intron   Promoter \n",
      "       270        105       8214         36        645       6915       1102 \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "GR_lead_taba  = makeGRangesFromDataFrame(data.frame(chr = LEDE_taba$Chr,start= LEDE_taba$Pos_Left,end = LEDE_taba$Pos_Right))\n",
    "GR_dep_taba  = makeGRangesFromDataFrame(data.frame(chr = LEDE_taba$dep_Chr,start= LEDE_taba$dep_Pos_Left,end = LEDE_taba$dep_Pos_Right))\n",
    "\n",
    "\n",
    "\n",
    "## get txdb of hg19\n",
    "txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene\n",
    "ah <- AnnotationHub()\n",
    "\n",
    "\n",
    "# annotate lead\n",
    "peakAnno <- annotatePeak(GR_lead_taba,tssRegion=c(-500, 500),\n",
    "                         TxDb=txdb, annoDb=\"org.Hs.eg.db\")\n",
    "\n",
    "anno_df = data.frame(peakAnno)\n",
    "anno_df$anno_short = sapply(anno_df$annotation,function(x){strsplit(x,split=\" \")[[1]][1]})\n",
    "\n",
    "# annotate dependents \n",
    "DAnno <- annotatePeak(GR_dep_taba,tssRegion=c(-500, 500),\n",
    "                         TxDb=txdb, annoDb=\"org.Hs.eg.db\")\n",
    "anno_dfd = data.frame(DAnno)\n",
    "anno_dfd$anno_short = sapply(anno_dfd$annotation,function(x){strsplit(x,split=\" \")[[1]][1]})\n",
    "\n",
    "head(anno_dfd)\n",
    "\n",
    "# add annotation to data frame\n",
    "LEDE_taba$lead_anno_tss500 = anno_df$anno_short\n",
    "LEDE_taba$dep_anno_tss500 = anno_dfd$anno_short\n",
    "print(table(LEDE_taba$lead_anno_tss500[!duplicated(LEDE_taba$Peak)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "8c5c5e8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Chr.1 Pos_Left Pos_Right                  pid  tad_id tad_chr tad_start\n",
      "1  chr1  1186158   1186820 chr1:1186158:1186820 tad_586    chr1   1030000\n",
      "2  chr1  1293775   1295119 chr1:1293775:1295119 tad_587    chr1   1255000\n",
      "3  chr1  1353131   1354479 chr1:1353131:1354479 tad_587    chr1   1255000\n",
      "4  chr1  1355362   1355777 chr1:1355362:1355777 tad_587    chr1   1255000\n",
      "5  chr1  1830444   1831294 chr1:1830444:1831294 tad_588    chr1   1710000\n",
      "6  chr1  2043563   2044398 chr1:2043563:2044398 tad_589    chr1   1860000\n",
      "  tad_end ab_id ab_chr ab_start  ab_end\n",
      "1 1235000    A1   chr1   700000 1900000\n",
      "2 1450000    A1   chr1   700000 1900000\n",
      "3 1450000    A1   chr1   700000 1900000\n",
      "4 1450000    A1   chr1   700000 1900000\n",
      "5 1840000    A1   chr1   700000 1900000\n",
      "6 2055000    B1   chr1  1900000 2100000\n",
      "      Peak Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7        7   1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222    213   1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.1  213   1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.2  213   1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287    277   1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296    286   1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "           Beta       SE  Log10_BF Pos_Left Pos_Right      GC Peak_Height\n",
      "7      0.886884 0.188153  3.085214   565005    565787 43.0396    13322207\n",
      "222    0.226942 0.024615 11.116547  2507916   2508523 59.2105        7712\n",
      "222.1  0.226942 0.024615 11.116547  2507916   2508523 59.2105        7712\n",
      "222.2  0.226942 0.024615 11.116547  2507916   2508523 59.2105        7712\n",
      "287   -0.654595 0.049849 19.692796  3819984   3820189 58.7379        1990\n",
      "296   -0.666789 0.041505 25.445000  4327332   4328105 53.8760        9809\n",
      "       P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7     0.985192 0.6912524  0.624720158        782      565396\n",
      "222   1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.1 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.2 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287   1.000000 0.9377345  0.020502304        205     3820086\n",
      "296   1.000000 0.9997858  0.004738315        773     4327718\n",
      "      Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr dep_Pos\n",
      "7                         110      0.1406650        8       1  569204\n",
      "222                       123      0.2026359      209       1 2508343\n",
      "222.1                     123      0.2026359      212       1 2508343\n",
      "222.2                     123      0.2026359      214       1 2508343\n",
      "287                        81      0.3951220      279       1 3820005\n",
      "296                        83      0.1073739      284       1 4327635\n",
      "         dep_RsID dep_Ref dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead  dep_Beta\n",
      "7     rs112660509       T       C 0.40125               1  0.4776132  0.251458\n",
      "222   rs200553981       G       C 0.35000               2  0.6725081 -0.134441\n",
      "222.1 rs200553981       G       C 0.35000               2  0.9016126 -0.205952\n",
      "222.2 rs200553981       G       C 0.35000               2  0.9806876  0.043439\n",
      "287    rs10909828       G       A 0.34500               2  0.9403693 -0.284635\n",
      "296     rs2101576       C       T 0.55000               2  0.9996411 -0.342060\n",
      "        dep_SE dep_Log10_BF dep_Pos_Left dep_Pos_Right  dep_GC dep_Peak_Height\n",
      "7     0.140264   -0.6716076       569112        569721 44.0984        29920971\n",
      "222   0.032999    1.3310342      2495782       2496858 64.3454             402\n",
      "222.1 0.046273    2.0625997      2499680       2500108 66.4336             660\n",
      "222.2 0.013210   -0.1630566      2508970       2510224 62.0717            6161\n",
      "287   0.051033    4.0437607      3826501       3826784 49.6479             556\n",
      "296   0.049869    6.4986698      4233102       4233752 52.8418             875\n",
      "      dep_P_caQTL      dep_PMR dep_P_eQTL_Coloc dep_Peak_width dep_Peak_center\n",
      "7        0.015921 6.650092e-07      0.015463689            609          569416\n",
      "222      0.032294 0.000000e+00      0.014451237           1076         2496320\n",
      "222.1    0.091908 0.000000e+00      0.042296138            428         2499894\n",
      "222.2    0.012672 0.000000e+00      0.001904334           1254         2509597\n",
      "287      0.997449 3.472557e-03      0.172399409            283         3826643\n",
      "296      0.970727 0.000000e+00      0.006162762            650         4233427\n",
      "      dep_Dist_SNP_to_peak_center dep_SNP_centrality LD_peak_dist\n",
      "7                             212          0.3481117         4020\n",
      "222                         12023         11.1737918        11900\n",
      "222.1                        8449         19.7406542         8326\n",
      "222.2                        1254          1.0000000         1377\n",
      "287                          6638         23.4558304         6557\n",
      "296                         94208        144.9353846        94291\n",
      "      enh_subcategory lead_anno_tss500 dep_anno_tss500 LEAD_TAD_id DEP_TAD_id\n",
      "7                  CM           Distal          Distal        <NA>       <NA>\n",
      "222                CM           Distal          Distal        <NA>       <NA>\n",
      "222.1              CM           Distal          Distal        <NA>       <NA>\n",
      "222.2              CM           Distal          Distal        <NA>       <NA>\n",
      "287                CM           Intron          Intron        <NA>       <NA>\n",
      "296                CM           Distal          Distal     tad_609    tad_609\n",
      "      LEAD_AB DEP_AB\n",
      "7        <NA>   <NA>\n",
      "222      <NA>   <NA>"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/site-packages/rpy2/robjects/pandas2ri.py:55: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.\n",
      "  for name, values in obj.iteritems():\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "222.1    <NA>   <NA>\n",
      "222.2    <NA>   <NA>\n",
      "287      <NA>   <NA>\n",
      "296        B1     B1\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "## add TAD information \n",
    "tadl = read.delim(\"/home/.../lead_peaks_GM_full_TAD_and_AB_overlap.bed\",header=T)\n",
    "tadd =read.delim(\"/home/.../dependent_peaks_GM_full_TAD_and_AB_overlap.bed\",header=T)\n",
    "print(head(tadl))\n",
    "lpid =sapply(1:length(LEDE_taba[,1]),function(x){paste(LEDE_taba[x,c(\"Chr\",\"Pos_Left\",\"Pos_Right\")],collapse=\":\")})\n",
    "dpid =sapply(1:length(LEDE_taba[,1]),function(x){paste(LEDE_taba[x,c(\"dep_Chr\",\"dep_Pos_Left\",\"dep_Pos_Right\")],collapse=\":\")})\n",
    "LEDE_taba$LEAD_TAD_id = tadl$tad_id[match(lpid,tadl$pid)]\n",
    "LEDE_taba$DEP_TAD_id = tadd$dep_tad_id[match(dpid,tadd$dep_pid)]\n",
    "LEDE_taba$LEAD_AB = tadl$ab_id[match(lpid,tadl$pid)]\n",
    "LEDE_taba$DEP_AB= tadd$dep_ab_id[match(dpid,tadd$dep_pid)]\n",
    "\n",
    "head(LEDE_taba)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d12bea0c",
   "metadata": {},
   "source": [
    "### add GC content of 700 bp used for scoring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "002cc40a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "#print(head(Gaf_LEDE_IPtab))\n",
    "lead_gr = makeGRangesFromDataFrame(data.frame(chr=LEDE_taba$Chr,start=LEDE_taba$Peak_center-350,\n",
    "                                              end=LEDE_taba$Peak_center+350))\n",
    "dep_gr = makeGRangesFromDataFrame(data.frame(chr=LEDE_taba$dep_Chr,start=LEDE_taba$dep_Peak_center-350,\n",
    "                                              end=LEDE_taba$dep_Peak_center+350))\n",
    "\n",
    "\n",
    "lead_seqs = as.character(getSeq(Hsapiens,lead_gr))\n",
    "#print(head(lead_seqs))\n",
    "dep_seqs = as.character(getSeq(Hsapiens,dep_gr))\n",
    "\n",
    "get_GC_cont=function(seq){\n",
    "    t=table(strsplit(seq,\"\")[[1]])[c('C','G')]\n",
    "    if (is.na(t[1]) & is.na(t[2])){\n",
    "        cg=0\n",
    "    }\n",
    "    if (!is.na(t[1]) | !is.na(t[2])){\n",
    "        cg = sum(t,na.rm=T)/nchar(seq)\n",
    "    }\n",
    "    return(cg)\n",
    "}\n",
    "\n",
    "LEDE_taba$caQTL_enh_GC_700bp = sapply(lead_seqs,get_GC_cont)\n",
    "LEDE_taba$dep_GC_700bp = sapply(dep_seqs,get_GC_cont)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4c192d9",
   "metadata": {},
   "source": [
    "### Add caQTL and no-effect label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "d89e78eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    caQTL no-effect \n",
      "     9676      7611 \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "LEDE_taba$enh_status_label = 'no-effect'\n",
    "LEDE_taba$enh_status_label[LEDE_taba$enh_subcategory =='CM' | LEDE_taba$enh_subcategory =='LOCAL'] = 'caQTL'\n",
    "print(table(LEDE_taba$enh_status_label[!duplicated(LEDE_taba$Peak)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49f70b37",
   "metadata": {},
   "source": [
    "### Add info about GM12878 SNPS and Allele-specific ATAC-seq counts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "a39f48d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  X.CHROM    POS         ID REF ALT GENOTYPE ATACseq_ref_count\n",
      "1       1 837214 rs72631888   G   C      0|1                 7\n",
      "2       1 851390 rs72631889   G   T      1|0                 5\n",
      "3       1 863511 rs28626846   T   C      1|0                 7\n",
      "4       1 875770  rs4970379   A   G      1|0                40\n",
      "5       1 876200 rs71628924   G   C      1|0                28\n",
      "6       1 900298 rs71628926   C   G      1|0                 8\n",
      "  ATACseq_alt_count  ATACseq_p\n",
      "1                 8 1.00000000\n",
      "2                 8 0.58105469\n",
      "3                 5 0.77441406\n",
      "4                33 0.48278856\n",
      "5                14 0.04355852\n",
      "6                 7 1.00000000\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "path_to_encode2 = '/data/.../'\n",
    "load(\"/home/.../Summary_table_SNPS_REF_ALT_counts_TFinHOCO_2021_04_08.RData\")\n",
    "\n",
    "file_list = list.files(path_to_encode2,full.names = TRUE)\n",
    "\n",
    "### load ATAC - counts at SNPs of GM12878 phased\n",
    "atac_tabn = read.table(paste0(path_to_encode2,'Combined_ATAC_R1R2_GM12878.ASB.summary.txt'),header=F)\n",
    "colnames(atac_tabn) = c(colnames(full_matf)[1:6],\"ATACseq_ref_count\",\"ATACseq_alt_count\",\"ATACseq_p\")\n",
    "print(head(atac_tabn))\n",
    "\n",
    "## GRanges of SNP positions in ATAC phased genome file\n",
    "GR_GM_SNPs = makeGRangesFromDataFrame(data.frame(chr = paste0(\"chr\",atac_tabn$X.CHROM),start = atac_tabn$POS,end=atac_tabn$POS))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "51c8168f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R \n",
    "## add the id of the lead RSiD in atac_tab\n",
    "LEDE_taba$LEAD_SNP_GMindex = match(LEDE_taba$RsID,atac_tabn$ID) \n",
    "## add GM info from full_matf as reference for double-checking\n",
    "LEDE_tabat = get_leadSNPs_in_GM_info(LEDE_taba,atac_tabn)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "7dab0ca1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Peak  Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7        7 chr1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222    213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.1  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.2  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287    277 chr1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296    286 chr1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "           Beta       SE  Log10_BF Chr Pos_Left Pos_Right Peak      GC\n",
      "7      0.886884 0.188153  3.085214   1   565005    565787    7 43.0396\n",
      "222    0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.1  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.2  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "287   -0.654595 0.049849 19.692796   1  3819984   3820189  277 58.7379\n",
      "296   -0.666789 0.041505 25.445000   1  4327332   4328105  286 53.8760\n",
      "      Peak_Height  P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7        13322207 0.985192 0.6912524  0.624720158        782      565396\n",
      "222          7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.1        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.2        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287          1990 1.000000 0.9377345  0.020502304        205     3820086\n",
      "296          9809 1.000000 0.9997858  0.004738315        773     4327718\n",
      "      Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr dep_Pos\n",
      "7                         110      0.1406650        8    chr1  569204\n",
      "222                       123      0.2026359      209    chr1 2508343\n",
      "222.1                     123      0.2026359      212    chr1 2508343\n",
      "222.2                     123      0.2026359      214    chr1 2508343\n",
      "287                        81      0.3951220      279    chr1 3820005\n",
      "296                        83      0.1073739      284    chr1 4327635\n",
      "         dep_RsID dep_Ref dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead  dep_Beta\n",
      "7     rs112660509       T       C 0.40125               1  0.4776132  0.251458\n",
      "222   rs200553981       G       C 0.35000               2  0.6725081 -0.134441\n",
      "222.1 rs200553981       G       C 0.35000               2  0.9016126 -0.205952\n",
      "222.2 rs200553981       G       C 0.35000               2  0.9806876  0.043439\n",
      "287    rs10909828       G       A 0.34500               2  0.9403693 -0.284635\n",
      "296     rs2101576       C       T 0.55000               2  0.9996411 -0.342060\n",
      "        dep_SE dep_Log10_BF dep_Chr dep_Pos_Left dep_Pos_Right dep_Peak  dep_GC\n",
      "7     0.140264   -0.6716076       1       569112        569721        8 44.0984\n",
      "222   0.032999    1.3310342       1      2495782       2496858      209 64.3454\n",
      "222.1 0.046273    2.0625997       1      2499680       2500108      212 66.4336\n",
      "222.2 0.013210   -0.1630566       1      2508970       2510224      214 62.0717\n",
      "287   0.051033    4.0437607       1      3826501       3826784      279 49.6479\n",
      "296   0.049869    6.4986698       1      4233102       4233752      284 52.8418\n",
      "      dep_Peak_Height dep_P_caQTL      dep_PMR dep_P_eQTL_Coloc dep_Peak_width\n",
      "7            29920971    0.015921 6.650092e-07      0.015463689            609\n",
      "222               402    0.032294 0.000000e+00      0.014451237           1076\n",
      "222.1             660    0.091908 0.000000e+00      0.042296138            428\n",
      "222.2            6161    0.012672 0.000000e+00      0.001904334           1254\n",
      "287               556    0.997449 3.472557e-03      0.172399409            283\n",
      "296               875    0.970727 0.000000e+00      0.006162762            650\n",
      "      dep_Peak_center dep_Dist_SNP_to_peak_center dep_SNP_centrality\n",
      "7              569416                         212          0.3481117\n",
      "222           2496320                       12023         11.1737918\n",
      "222.1         2499894                        8449         19.7406542\n",
      "222.2         2509597                        1254          1.0000000\n",
      "287           3826643                        6638         23.4558304\n",
      "296           4233427                       94208        144.9353846\n",
      "      LD_peak_dist enh_subcategory lead_anno_tss500 dep_anno_tss500 LEAD_TAD_id\n",
      "7             4020              CM           Distal          Distal        <NA>\n",
      "222          11900              CM           Distal          Distal        <NA>\n",
      "222.1         8326              CM           Distal          Distal        <NA>\n",
      "222.2         1377              CM           Distal          Distal        <NA>\n",
      "287           6557              CM           Intron          Intron        <NA>\n",
      "296          94291              CM           Distal          Distal     tad_609\n",
      "      DEP_TAD_id LEAD_AB DEP_AB caQTL_enh_GC_700bp dep_GC_700bp causal_label\n",
      "7           <NA>    <NA>   <NA>          0.4136947    0.4407989    no-effect\n",
      "222         <NA>    <NA>   <NA>          0.5962910    0.6447932    no-effect\n",
      "222.1       <NA>    <NA>   <NA>          0.5962910    0.6604850    no-effect\n",
      "222.2       <NA>    <NA>   <NA>          0.5962910    0.6276748    no-effect\n",
      "287         <NA>    <NA>   <NA>          0.5734665    0.5492154    no-effect\n",
      "296      tad_609      B1     B1          0.5335235    0.5235378    no-effect\n",
      "      enh_status_label LEAD_SNP_GMindex     GM_RsID GM_REF GM_ALT  GM_Pos\n",
      "7                caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "222              caQTL              737 rs200553981      G      C 2508343\n",
      "222.1            caQTL              737 rs200553981      G      C 2508343\n",
      "222.2            caQTL              737 rs200553981      G      C 2508343\n",
      "287              caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "296              caQTL             2125   rs2101576      C      T 4327635\n",
      "      GM_phasing Dist_GM_SNP_to_peak_center\n",
      "7           <NA>                         NA\n",
      "222          0|1                        123\n",
      "222.1        0|1                        123\n",
      "222.2        0|1                        123\n",
      "287         <NA>                         NA\n",
      "296          1|0                         83\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "print(head(LEDE_tabat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "f67cf0bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       \n",
      "           CM LOCAL no-effect\n",
      "  FALSE  2649 10848     10960\n",
      "  TRUE   1259  5223      5180\n",
      "      Peak  Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7        7 chr1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222    213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.1  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.2  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287    277 chr1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296    286 chr1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "           Beta       SE  Log10_BF Chr Pos_Left Pos_Right Peak      GC\n",
      "7      0.886884 0.188153  3.085214   1   565005    565787    7 43.0396\n",
      "222    0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.1  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.2  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "287   -0.654595 0.049849 19.692796   1  3819984   3820189  277 58.7379\n",
      "296   -0.666789 0.041505 25.445000   1  4327332   4328105  286 53.8760\n",
      "      Peak_Height  P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7        13322207 0.985192 0.6912524  0.624720158        782      565396\n",
      "222          7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.1        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.2        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287          1990 1.000000 0.9377345  0.020502304        205     3820086\n",
      "296          9809 1.000000 0.9997858  0.004738315        773     4327718\n",
      "      Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr dep_Pos\n",
      "7                         110      0.1406650        8    chr1  569204\n",
      "222                       123      0.2026359      209    chr1 2508343\n",
      "222.1                     123      0.2026359      212    chr1 2508343\n",
      "222.2                     123      0.2026359      214    chr1 2508343\n",
      "287                        81      0.3951220      279    chr1 3820005\n",
      "296                        83      0.1073739      284    chr1 4327635\n",
      "         dep_RsID dep_Ref dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead  dep_Beta\n",
      "7     rs112660509       T       C 0.40125               1  0.4776132  0.251458\n",
      "222   rs200553981       G       C 0.35000               2  0.6725081 -0.134441\n",
      "222.1 rs200553981       G       C 0.35000               2  0.9016126 -0.205952\n",
      "222.2 rs200553981       G       C 0.35000               2  0.9806876  0.043439\n",
      "287    rs10909828       G       A 0.34500               2  0.9403693 -0.284635\n",
      "296     rs2101576       C       T 0.55000               2  0.9996411 -0.342060\n",
      "        dep_SE dep_Log10_BF dep_Chr dep_Pos_Left dep_Pos_Right dep_Peak  dep_GC\n",
      "7     0.140264   -0.6716076       1       569112        569721        8 44.0984\n",
      "222   0.032999    1.3310342       1      2495782       2496858      209 64.3454\n",
      "222.1 0.046273    2.0625997       1      2499680       2500108      212 66.4336\n",
      "222.2 0.013210   -0.1630566       1      2508970       2510224      214 62.0717\n",
      "287   0.051033    4.0437607       1      3826501       3826784      279 49.6479\n",
      "296   0.049869    6.4986698       1      4233102       4233752      284 52.8418\n",
      "      dep_Peak_Height dep_P_caQTL      dep_PMR dep_P_eQTL_Coloc dep_Peak_width\n",
      "7            29920971    0.015921 6.650092e-07      0.015463689            609\n",
      "222               402    0.032294 0.000000e+00      0.014451237           1076\n",
      "222.1             660    0.091908 0.000000e+00      0.042296138            428\n",
      "222.2            6161    0.012672 0.000000e+00      0.001904334           1254\n",
      "287               556    0.997449 3.472557e-03      0.172399409            283\n",
      "296               875    0.970727 0.000000e+00      0.006162762            650\n",
      "      dep_Peak_center dep_Dist_SNP_to_peak_center dep_SNP_centrality\n",
      "7              569416                         212          0.3481117\n",
      "222           2496320                       12023         11.1737918\n",
      "222.1         2499894                        8449         19.7406542\n",
      "222.2         2509597                        1254          1.0000000\n",
      "287           3826643                        6638         23.4558304\n",
      "296           4233427                       94208        144.9353846\n",
      "      LD_peak_dist enh_subcategory lead_anno_tss500 dep_anno_tss500 LEAD_TAD_id\n",
      "7             4020              CM           Distal          Distal        <NA>\n",
      "222          11900              CM           Distal          Distal        <NA>\n",
      "222.1         8326              CM           Distal          Distal        <NA>\n",
      "222.2         1377              CM           Distal          Distal        <NA>\n",
      "287           6557              CM           Intron          Intron        <NA>\n",
      "296          94291              CM           Distal          Distal     tad_609\n",
      "      DEP_TAD_id LEAD_AB DEP_AB caQTL_enh_GC_700bp dep_GC_700bp causal_label\n",
      "7           <NA>    <NA>   <NA>          0.4136947    0.4407989    no-effect\n",
      "222         <NA>    <NA>   <NA>          0.5962910    0.6447932    no-effect\n",
      "222.1       <NA>    <NA>   <NA>          0.5962910    0.6604850    no-effect\n",
      "222.2       <NA>    <NA>   <NA>          0.5962910    0.6276748    no-effect\n",
      "287         <NA>    <NA>   <NA>          0.5734665    0.5492154    no-effect\n",
      "296      tad_609      B1     B1          0.5335235    0.5235378    no-effect\n",
      "      enh_status_label LEAD_SNP_GMindex     GM_RsID GM_REF GM_ALT  GM_Pos\n",
      "7                caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "222              caQTL              737 rs200553981      G      C 2508343\n",
      "222.1            caQTL              737 rs200553981      G      C 2508343\n",
      "222.2            caQTL              737 rs200553981      G      C 2508343\n",
      "287              caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "296              caQTL             2125   rs2101576      C      T 4327635\n",
      "      GM_phasing Dist_GM_SNP_to_peak_center lead_peak_id dep_peakid\n",
      "7           <NA>                         NA            7          8\n",
      "222          0|1                        123          213        209\n",
      "222.1        0|1                        123          213        212\n",
      "222.2        0|1                        123          213        214\n",
      "287         <NA>                         NA          277        279\n",
      "296          1|0                         83          286        284\n",
      "      depGMSNP_RsID depGMSNP_REF depGMSNP_ALT depGMSNP_Pos depGMSNP_phasing\n",
      "7              <NA>         <NA>         <NA>           NA             <NA>\n",
      "222      rs10910090            C            T      2496214              1|0\n",
      "222.1     rs4385650            G            C      2499898              1|0\n",
      "222.2          <NA>         <NA>         <NA>           NA             <NA>\n",
      "287            <NA>         <NA>         <NA>           NA             <NA>\n",
      "296            <NA>         <NA>         <NA>           NA             <NA>\n",
      "      depGMSNP_dist_to_dep_peak_center depGMSNP_refcount depGMSNP_altcount\n",
      "7                                   NA                NA                NA\n",
      "222                                106                12                13\n",
      "222.1                                4                16                33\n",
      "222.2                               NA                NA                NA\n",
      "287                                 NA                NA                NA\n",
      "296                                 NA                NA                NA\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAMAAABKCk6nAAAC9FBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpcXFxeXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+RkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqcnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///8v2GqzAAAWkklEQVR4nO2deXwUZZrHeyFg1AABwiEGCKyEw4wcsoAQEQVBMI4oIHiMooMiyzh4DK7Kgo7HKIo3IiA4rC6gK66uOiKDShyJgsMikqCJMiwKYg7EmDt5/tnqrupUU8dr6u2q7vd9+X0/TCdV9fT79vTX7v5Vut+nQwSUJpTsGwCCBYIVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVJw7BP2wAAvBKXVCC11+1AiSfMcWBCX6G/7rAN66HYLXxRfCxJoedECwEcQuueiD7pFBKv3trrQcgWAjiFnzd1Pzy+vK/zbzeegCChSBuwRlVkR8NvawHIFgI4hY89LXIj78Osx6AYCGIW/D2nr+aecOsIZk7rQcgWAjiT9ENm1c+/PyWBtt+CBYCf86Dyyrt+yBYCOIWnFdNxSNatZl82HoAgoUgbsGhShq/sKp20XTrAQ7Bb//JmZvvdTnwrecpHKlb6jL+bJf9q/2Zl3a6jH/7HS4H/u55Cj8Ep9drr8Rdm/dsXRhhvO3M+Bcp2uzMZStcDhzzPIUjTX91GX+wy/5P/JmXDrqMf9sfXA54/y86fsElNKyEqLBf854fdkS48nLPN8aNW3b5NpQ3xiVp3lVrfRsqbsFjM9PSp9C2rk9aD/x+Jv+tsnDL//o2lDeSJXj1i74N5UOKri4soIIttt0+Ct5v+0N3gvgqSfOWlvk2lE9vF9bbd/koGPATt+DD149bsbNf63P3Ww/4KPhvP/o2lDfeSdK8X7KteCFuwRdfvuH8Hu9V3JdnPeDnazBCFjdxC06voDeHEDV2tR6AYH5EEpxZSKVawjrcz3oAKZofkVL0M91u0y5XDHrYegApmh+hUvSe17WLp9+y7UeKFoLgPlWJFM2PSCnaFYQsfkQKWa5AMD8nnGCkaG6kEIwUzY8UggE/UghGiuZHCsEIWfxAMAsIZoAUzQ9SdIJAimaAFC0EUghGiuZHCsEIWfxAMAsIZoAUzQ9SdIJAimaAFC0EUghGiuZHCsEIWfxAMAsIZoAUzQ9SdIJAimaAFC0EUghGiuYnuG6zCFn8iBSyXLvNQjA/Igl27TaLFM2PSCnatdssUjQ/IqVo126zSNFCEFy3WaRofoRK0W7dZhGy+BEpZIXBaZLPiCQ4IadJSNHcSHGahBTNTwCnSa9OiNBrfHy3DPhCAKdJDeUR5iJFcyNSik7EaRJCFj84TWIBwVGC7TaLFM2PFN1mkaL5kaLbLOBHim6zSNH8yNFtFiGLGym6zUIwP1J0m0WK5keKT1UiRfMjhWDAjxSCkaL5kUIwQhY/EMwCghkgRfODFJ0gkKIZIEULgRSCkaL5kUIwQhY/EMwCghkgRfODFJ0gkKIZIEULgRSCkaL5kUIwQhY/EMwCghkgRfODFJ0gkKIZIEULgRSCkaL5kUIwQhY/EMwCghkgRfODFJ0gkKJ1gm6jBPiRoo0SUjQ/UrRRQsjiR4pusxDMjxTdZpGi+ZGijRJSND9StFEC/PhzHlxWad+HFM2PSCk6r5qKR7RqM/mw9QBCFj8ihaxQJY1fWFW7aLr1AATzI5jg9Hrtldhso/SK0W32grhvXBSkaH7iF1xCw0qICoNso4QUzU/cgsdmpqVPoW1dn7QeQIoWgpYK/ia/2OkdhTDVhQVUsMW2Gyman4Sn6H1ndxqcMbjIuaYicnnUuhshi5+Eh6wx99VTwwOjnSp2D2rVdxNRje3JHIL5Sbjg9uF+0A0ZThWjl9Zt7Z7vTXDeIG907+etvs96D/cAixMmRZ/9qnbx2iVOFR0aiV4fWOdJ8LjPg+Wxpz3cAyxOmBQ9LTRixojQpKuv/thWcda72sXMeWoKVoCWCV4XxdbXnd46NbeMKkacraTgEyZFEzWVupX838uVRHUbF1r3qyD4hAlZ/zjv5JyPzrc/fFlAMD8JFzzj7uqc+rsv9DSyiIKPlnsj12P9z57uIXcSnqK7NVAO1XT0NLKAgnf1GOWNIR7rOzd6uotcSXiKPitfE/zRIE8jCyi4YHbAE4+0f7Al2bRM8NYuV2VM7/KGp5Eh2KR4szfWrPV4hRLX+7qFKfqHVUuWH2BW2oBgk+tmL/CE1/prf+t6XydlbdKJJ/idYOd9Mz7BGTpdspilViDYRHDBR48+NDH/0M5LlzNLrUCwieCCiXr/pF1U9WGWWoFgE+EFn75XuyjqySy1AsEmwgte1n3xmiWnLWOWWoFgE+EF09Z50+a/z6y0AcEm4gs2uI1ZfDwQbCKN4JuYxccDwSYQ7C8QbALBfgLBOhDsMxCsA8E2INgT8greR7Ub1jZQIbP4eCDYRHjB97SnxUOG38wstQLBJsIL7vxNY+eism7MUisQbCK84E4VH/Sn6k7MUisQbCK84DlDs577btylzFIrEGwivOCGDRsaDyz9ya3Ka7dZCPaZeAV/8bL7Al2ebrMQ7DNxCl6ROjj1KbcKnm6zEOwzcQrOeoU2dXer4Ok2C8E+E6fg1jVUl+JWwdNtFoJ9Jk7BYbmpbhU83WYh2GfifQTv2LGjrfa/HY4lHN1mIdhn/Pnge4ZjE5YwOE0ykFQwE5wmxaCiYPtpUoO+GHouBJt3ktiCJ0RxqrCfJkW7zY53HRWCfSZOwZsiPHDqEKcKnCbFIKngMFV3dXzM8dbjNCkGeQW/3fdSt/XfOE0ykVXwdzN7bnKv8d5tFoJ9Jk7Bz3S+1eFbVQx4us1CsM/EKTjUOk3HqYKn2ywE+0ycgg9HcaqQqdssBNswtP1E9MkLzp0MZeo2C8E2ItqKz5xGH7bOzXjPqUKmbrMQbCOi7Yr7G+ncx2h9rmOJRN1mIdiG/omOWqo46SiVpTNLrUCwieCCO1XQxnO0M13Xj+04AsEmggueuOzYuQ8SPXsxs9QKBJsILvjzXv+UU0mTO37JLLUCwSaCC6bGg9rFGsfTYHcg2ER0wc2gyw4b6QVjATgbCNaBYJ+BYB0ItgHBnoBgHQj2GQjWgWAbEcH7o6DLzi8gqeBQ23QdZqkVCDYRXPDNvSY+f4RZ5gQEmwgumJryF/Qd9/S3zEobEGwiuuAwn96ZPYZZagWCTSQQ3Pj+vJ7Sf/soBNvQBTdsmZt50UrXr4h2BoJNBBc85/QpL5Qzy5yAYBPBBYfaMD747goEmwgu+GAUZqkVCDYRXDDV76X7lyxZspJZagWCTQQXfCjnRkr9w7CeHzJLrUCwieCCZ88J98lqvHY1s9QKBJsILrjHZ5FGaAf+xa0KbZQMJBV8ys/hL22gqg5OFWijFIOkgodsifz+P6Mcbzu6zZpIKvjF/l9plx/1WudUgW6zMUgqmJa3y52ak7bUsQJtlGKQVTAd2vT4Jpd3hNFGKQZpBUdxWtmANkom0gt2+9BdmUMXHgg2kVlwXjUVj2jVZrJtaRoEm8gsOFRJ4xdW1S6a3rznhx0RrrzcdVQI9pmABafXa6/EXZv3bF0YYfgk11Fz3w6Wu1y+IwaC7fqO33QSXELDSogK+1kPMJ6izwgFzDLneSHYrk//UWtcOggem5mWPoW2dX3SegCCTYQXnBO+qMl0XtlQXVhABVtsuxmC8RTtM3EKnpISSgkzlVlqBSnaRHDBRJcwi5yBYBPhBUe/4t2Bb4sMrAcg2ER4wYyveH+odUr3CNYDIgo+c2KwZMkqmPUV7xtznK8qouD0gNN7hqyCWV/xXnm181VFFDxidrD0k1UwvuK9ZUj7GvwLX/HuCASbXPNCsCf+K2e73tdxf7WdKxBsMjng1/7Qr13v65YIHl4wXIdZagWCTQQXvP1H/Q1A5+8PdgOCTfAU7S/CCRY8ZBlfD90li1lqBYJNBBd89OhDE/MP7bx0ObPUCgSbCC6YqHf4DKmqD7PUCgSbCC/49L3aRVFPZqkVCDYRXvCy7ovXLDnN5XMSLkCwifCCaeu8afPfZ1bagGAT8QVzAMEmwgvenJsThllqBYJNhBfc649fOH1ogwkEmwgvuItTj4ZfAIJNhBd85XLvhiHYRHjBY0Kd+g8YMIBZagWCTYQXvEeHWWoFgk2EF8wDBJsILhhv+LcUSQXjDf+WIqlgjTL9a3W8AMEmwgv+Xahnlgaz1AoEmwgvuN1uZpUjEGwivOCRZcwqRyDYRHjBf8m85X4Ntyp0mzWQVvA5o//tHg2nCnSbjUFawR2rXSvQbTYGaQVfsMn1zQZ0m41BWsFDQ90HuLzZgG6zMUgrmPFmA7rNxiCtYBboNmuipGDCaVIzKgrGaVIMKgq2nybt2xDhYvc1qxDsM0EKtp8m7V4RYXye63Ug2GeCFIzTpBhUFIzTpBhUFIzTpBiUFOwGBJtAsL9AsEncgtFtNgYVBUvVbRaCrbTgKVqmbrMQbKUFgmXqNgvBVhCyPAHBOhDsMxCsA8E2INgTEKwDwT4DwToQbAOCPQHBOhDsMxCsA8E2INgTEKwDwT4DwToQbAOCPQHBOhDsMxCsA8E2INgTEKwDwT4DwToQbAOCPQHBOrdO8Ebvczxe4V3neSHYRjCCvXLLLn/GgWAbYgi+/Qt/xoFgG2IIrvdpHAi2IYZgv0ie4NkLAuVayQW/fMifcT45a7Y3LvBY38NF8N4N3vjjAx6vsM/1/7MUgv0KWfU7PHK2x3qOrsuOrFrr00AnmGDPjEvSvCeaYL9StGemJGneF1/ybSgpBPuVoqWZuLHRt6GS0ukOJI6kdLrzil8p2jOPJ2nej7b7NlQAne4MELL4ESlk8TQE9woE85OUTndeQYrmJymd7ryCFM1PAJ3uPlgYYfikOG8a8AN/zoPLKmM2juh/trvycv5bZQEpmp+4BedVU/GIVm0mH7YeQMjiR6SQFaqk8QurahdNtx6AYH4EE5yuZZGGrtYDSNH8iJSiQyU0rISosJ/1AFI0PyKl6LGZaelTaFvXJ60H8LdoIfAhRVcXFlDBFttu/C2aH5FStCsIWfyIFLJcgWB+TjTBSNH8SCEYKZofKQQDfqQQjBTNjxSCEbL4gWAWEMwAKZofpOgEgRTNAClaCKQQjBTNjxSCEbL4gWAWEMwAKZofpOgEgRTNAClaCKQQjBTNjxSCEbL4gWAWEMwAKZofpOgEgRTNAClaCKQQjBTNjxSCEbL4gWAWEMwAKZofpOgEgRTNAClaCKQQjBTNjxSCEbL4gWAWEMwAKZofpOgEgRTNAClaCKRoCI4UzY8UDcERsvhBQ3AWEOzUEPz1GRGyJ3q+MctnONN/ksuBIs9TOFIzy2X801z23+bPvPS2y/gjRroc+IvnKQJoCF5XHmHlE55vTHW5M0dc9pc7vfbz8KPHiSt/ecgWUe8yfmmpywHvsT64huDrn/F8Y4D/BNAQ3ACChSC482AIFgIIVhwIVhwIVhwIVhwIVhwIVhwIVhwIVpzgBL8zeIJf9EjvlBxSkzRvuwG+3Xf9vw1KsI/g7cLAgOCkAMEJAoITAwQHBgQnBQhOEBCcGBbsTtLEE5I075p1iZpJDME/+fXZK68cS9K8NbYPIQeFGIJBYECw4kCw4kCw4kCw4kCw4kCw4kCw4iRZ8Mohp/Z9mOhgaE546+bQ/gDnqgwd1X9ZN7r9wAWV5uwa/9ytIaYgISzuWar9O25XEB1Dkit4WU5BzSd9n6GDbU7X7uCmXm0TIXhx9pafi2b9qrZ5dqJPM7O2JFpwh7Lwv+NICWCapAou61CoXa6/gg6mTdHu4PyJ6QkQXNKuJPxj9LLm2YluvXPhb4MRvGfMVV2y3iPamJ0x7Qd91/uDu884QDNbDZzQamCpvkEvndHtxupLQjn+/wUzqYI3n2X8cjDthblEC55PhODVF0W2Vk1unp2aMvfu6lwXjODQs0335VJxx09q5lwV2VPaJb/hrvOI0iq1f8bG3m6F349ard4j+PlJxi8H08ozG6nv94kQ/O+Rl3va0r95dvpgBNGAt4IR3KmJvsihR68n+v6kyFsqq68kqj65SRdsbCy5leiLj9UT/E7kMVT5XN3BNLro/U/HUSIEr9KjzJq85tlpbmpGRtvfBCO4P1FRDt1+r/Z7aunKM87YsCQjR6NUF2xs3KA3S1BO8JG0vdrla720RzCtnHfn0wkR/GW7yCTnPdY8e32XIqLd7WsCETwgIvjRG7Tn5raR7mfLpxPVf2Q8RRsbdy8k2rFRPcH04KCC2u29Hw0LLu2Z/W1CBNNdg7bVfH3NqPrm2d/JDu/uvSlAwV922lV/8xWRPQe6fFx7z/mGYGNjZ4/i0tzHKKXG9/mTLLjpiYGn9FvaGBZMF46moAWflJqa+oL27Dyqfdb8MnP2a+8IH54/q7nARwzBtD67y2VH9F1vDEyfsN8QHN1Y2afztbU0NVOxFA2CB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IVB4IV50QTnKwvKk4aEgmOfugxo3WKxn5Kb51y8qitMTsMmpdwrbv6uOtr+4880net9Qug9uS4zmgZwDK8HMgo2Oialr6Hfvpzh8/MHQbNny62+Emhij6Lh112jWVYr4LdP7wcwIde40duwUR3zbIKDi/heiqr3cgiWjd1WqdztCJ9S9v/nxc35B29ILr8y1gXtmfQnb16vxu9dsONnTMWUX4u0fbc5gEiS8PMYWqjK8ay+t6TG11Xtmvc3Wdf9zjR3bcm7h5pCbIL/jDb4RFc0e7zmptupHWhV+vuPbPJ2NL2f5xZkkfNy7+MdWF7Qs81PTEmeuWNQ8u/ab8vKlgfQF8aZg5jXP/r7vvKhzavK9t16iN1r15IlJOfoLujhcgoOPyS28EQXJRm7jBIodrvqWL+LFo3kqi+41fGVvi59d6Mbo+WRZd/GevC9nRuot3NT9Prsz5spOZHsD6AvjTMHMa4/tJ5RCtyo+vKdnVspGPtKou7Nyb8jmEio+DYR/C2fg6P4Po7zhw7SRM8Q9sanG9sRV48fxwyLbvKWPFlrAvbc2bs63Djs9mdl9SGBX+cGx1AXxpmDmNc//cPEr2ZG11XtmuQ9uPC1x+9Kfg7whOyC14000HwS8NK6UVN8CjtNTXjH8aWtv+D/IY8Omu7seLLWFUSlmsKPvAdFQ7ZkD+a6I3c6AD60jBzGOP6D/4r0crc6LqyXeEhnrgxd3Oi7o4WIrfgqg3pDim65qkL6o7mXq69hP5344OjyNjS9q+4rCGvts8hY8WXo+BHxlWXD1+7++QDNRfmRgfQl4aZwxjX33va1z+OODe6riwiuLhLhvdveA4WqQTra8My2mo/UjdQetvUU0Zu1Q5EdxhMzTwyoVvuf3T/r3UXj+183td0VN/S9leMntiz/5Lo8i9Hwcd+3bHr3PqmOe3bzdUE6wPoS8PMYWqNFWOreuQsmRpdVxYRTAN+k/j7hY1Egn2g8e/nlPs22O7bm+h3Dxy/L/d134b3CXUEF94U4TZmUeOf47n28TTM7ztUOz2Oof6zztVu1clCHcEC8Npp65N9E2xAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOJAsOL8P3NoNMiAkKlfAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "# add GM RSID if in dependent peak \n",
    "dep_GM_snp_info_tab = get_SNPs_in_deppeak_info_plus(LEDE_tabat,atac_tabn)\n",
    "## above should be computing the sum of all columns in atac_tabn for the GM dep snp (partially redundant bc done again below..)\n",
    "\n",
    "# choose the SNP in the dep peak that is closest to the peak center\n",
    "mindist_GM_SNPs_dep = sapply(dep_GM_snp_info_tab,function(i){\n",
    "  if(is.na(i$depGMSNP_dist_to_dep_peak_center[1])){\n",
    "    return(NA)} \n",
    "  if (!is.na(i$depGMSNP_dist_to_dep_peak_center[1])){\n",
    "    return(min(i$depGMSNP_dist_to_dep_peak_center,na.rm=T))}\n",
    "})\n",
    "boxplot(mindist_GM_SNPs_dep~LEDE_tabat$enh_subcategory)\n",
    "print(table(!is.na(mindist_GM_SNPs_dep),LEDE_tabat$enh_subcategory))\n",
    "\n",
    "# add the mindist chosen dep GM SNP to the master table\n",
    "LEDE_taba_f= cbind(LEDE_tabat,do.call(rbind,lapply(dep_GM_snp_info_tab,function(x){\n",
    "  if(is.na(x$depGMSNP_dist_to_dep_peak_center[1])){\n",
    "    return(data.frame(x[1,]))}\n",
    "  if (!is.na(x$depGMSNP_dist_to_dep_peak_center[1])){\n",
    "    return(data.frame(x[which.min(x$depGMSNP_dist_to_dep_peak_center),]))}\n",
    "})))\n",
    "\n",
    "print(head(LEDE_taba_f))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "d6019266",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Peak  Chr     Pos        RsID Ref Alt     AF Inside_Peak    P_Lead\n",
      "7        7 chr1  565286   rs1578391   C   T 0.5978           1 0.9985614\n",
      "222    213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.1  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "222.2  213 chr1 2508343 rs200553981   G   C 0.3500           1 0.9999242\n",
      "287    277 chr1 3820005  rs10909828   G   A 0.3450           1 0.9999658\n",
      "296    286 chr1 4327635   rs2101576   C   T 0.5500           1 0.9997817\n",
      "           Beta       SE  Log10_BF Chr Pos_Left Pos_Right Peak      GC\n",
      "7      0.886884 0.188153  3.085214   1   565005    565787    7 43.0396\n",
      "222    0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.1  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "222.2  0.226942 0.024615 11.116547   1  2507916   2508523  213 59.2105\n",
      "287   -0.654595 0.049849 19.692796   1  3819984   3820189  277 58.7379\n",
      "296   -0.666789 0.041505 25.445000   1  4327332   4328105  286 53.8760\n",
      "      Peak_Height  P_caQTL       PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "7        13322207 0.985192 0.6912524  0.624720158        782      565396\n",
      "222          7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.1        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "222.2        7712 1.000000 0.9800680  0.053831988        607     2508220\n",
      "287          1990 1.000000 0.9377345  0.020502304        205     3820086\n",
      "296          9809 1.000000 0.9997858  0.004738315        773     4327718\n",
      "      Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr dep_Pos\n",
      "7                         110      0.1406650        8    chr1  569204\n",
      "222                       123      0.2026359      209    chr1 2508343\n",
      "222.1                     123      0.2026359      212    chr1 2508343\n",
      "222.2                     123      0.2026359      214    chr1 2508343\n",
      "287                        81      0.3951220      279    chr1 3820005\n",
      "296                        83      0.1073739      284    chr1 4327635\n",
      "         dep_RsID dep_Ref dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead  dep_Beta\n",
      "7     rs112660509       T       C 0.40125               1  0.4776132  0.251458\n",
      "222   rs200553981       G       C 0.35000               2  0.6725081 -0.134441\n",
      "222.1 rs200553981       G       C 0.35000               2  0.9016126 -0.205952\n",
      "222.2 rs200553981       G       C 0.35000               2  0.9806876  0.043439\n",
      "287    rs10909828       G       A 0.34500               2  0.9403693 -0.284635\n",
      "296     rs2101576       C       T 0.55000               2  0.9996411 -0.342060\n",
      "        dep_SE dep_Log10_BF dep_Chr dep_Pos_Left dep_Pos_Right dep_Peak  dep_GC\n",
      "7     0.140264   -0.6716076       1       569112        569721        8 44.0984\n",
      "222   0.032999    1.3310342       1      2495782       2496858      209 64.3454\n",
      "222.1 0.046273    2.0625997       1      2499680       2500108      212 66.4336\n",
      "222.2 0.013210   -0.1630566       1      2508970       2510224      214 62.0717\n",
      "287   0.051033    4.0437607       1      3826501       3826784      279 49.6479\n",
      "296   0.049869    6.4986698       1      4233102       4233752      284 52.8418\n",
      "      dep_Peak_Height dep_P_caQTL      dep_PMR dep_P_eQTL_Coloc dep_Peak_width\n",
      "7            29920971    0.015921 6.650092e-07      0.015463689            609\n",
      "222               402    0.032294 0.000000e+00      0.014451237           1076\n",
      "222.1             660    0.091908 0.000000e+00      0.042296138            428\n",
      "222.2            6161    0.012672 0.000000e+00      0.001904334           1254\n",
      "287               556    0.997449 3.472557e-03      0.172399409            283\n",
      "296               875    0.970727 0.000000e+00      0.006162762            650\n",
      "      dep_Peak_center dep_Dist_SNP_to_peak_center dep_SNP_centrality\n",
      "7              569416                         212          0.3481117\n",
      "222           2496320                       12023         11.1737918\n",
      "222.1         2499894                        8449         19.7406542\n",
      "222.2         2509597                        1254          1.0000000\n",
      "287           3826643                        6638         23.4558304\n",
      "296           4233427                       94208        144.9353846\n",
      "      LD_peak_dist enh_subcategory lead_anno_tss500 dep_anno_tss500 LEAD_TAD_id\n",
      "7             4020              CM           Distal          Distal        <NA>\n",
      "222          11900              CM           Distal          Distal        <NA>\n",
      "222.1         8326              CM           Distal          Distal        <NA>\n",
      "222.2         1377              CM           Distal          Distal        <NA>\n",
      "287           6557              CM           Intron          Intron        <NA>\n",
      "296          94291              CM           Distal          Distal     tad_609\n",
      "      DEP_TAD_id LEAD_AB DEP_AB caQTL_enh_GC_700bp dep_GC_700bp causal_label\n",
      "7           <NA>    <NA>   <NA>          0.4136947    0.4407989    no-effect\n",
      "222         <NA>    <NA>   <NA>          0.5962910    0.6447932    no-effect\n",
      "222.1       <NA>    <NA>   <NA>          0.5962910    0.6604850    no-effect\n",
      "222.2       <NA>    <NA>   <NA>          0.5962910    0.6276748    no-effect\n",
      "287         <NA>    <NA>   <NA>          0.5734665    0.5492154    no-effect\n",
      "296      tad_609      B1     B1          0.5335235    0.5235378    no-effect\n",
      "      enh_status_label LEAD_SNP_GMindex     GM_RsID GM_REF GM_ALT  GM_Pos\n",
      "7                caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "222              caQTL              737 rs200553981      G      C 2508343\n",
      "222.1            caQTL              737 rs200553981      G      C 2508343\n",
      "222.2            caQTL              737 rs200553981      G      C 2508343\n",
      "287              caQTL               NA        <NA>   <NA>   <NA>      NA\n",
      "296              caQTL             2125   rs2101576      C      T 4327635\n",
      "      GM_phasing Dist_GM_SNP_to_peak_center lead_peak_id dep_peakid\n",
      "7           <NA>                         NA            7          8\n",
      "222          0|1                        123          213        209\n",
      "222.1        0|1                        123          213        212\n",
      "222.2        0|1                        123          213        214\n",
      "287         <NA>                         NA          277        279\n",
      "296          1|0                         83          286        284\n",
      "      depGMSNP_RsID depGMSNP_REF depGMSNP_ALT depGMSNP_Pos depGMSNP_phasing\n",
      "7              <NA>         <NA>         <NA>           NA             <NA>\n",
      "222      rs10910090            C            T      2496214              1|0\n",
      "222.1     rs4385650            G            C      2499898              1|0\n",
      "222.2          <NA>         <NA>         <NA>           NA             <NA>\n",
      "287            <NA>         <NA>         <NA>           NA             <NA>\n",
      "296            <NA>         <NA>         <NA>           NA             <NA>\n",
      "      depGMSNP_dist_to_dep_peak_center depGMSNP_refcount depGMSNP_altcount\n",
      "7                                   NA                NA                NA\n",
      "222                                106                12                13\n",
      "222.1                                4                16                33\n",
      "222.2                               NA                NA                NA\n",
      "287                                 NA                NA                NA\n",
      "296                                 NA                NA                NA\n",
      "      ATAC_caQTL_ref ATAC_caQTL_alt ATAC_caQTL_more ATAC_caQTL_less\n",
      "7                 NA             NA              NA              NA\n",
      "222                6             10              10               6\n",
      "222.1              6             10              10               6\n",
      "222.2              6             10              10               6\n",
      "287               NA             NA              NA              NA\n",
      "296               69              1              69               1\n",
      "      ATAC_dep_ref ATAC_dep_alt ATAC_dep_more ATAC_dep_less\n",
      "7               NA           NA            NA            NA\n",
      "222             12           13            12            13\n",
      "222.1           16           33            16            33\n",
      "222.2           NA           NA            NA            NA\n",
      "287             NA           NA            NA            NA\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "296             NA           NA            NA            NA\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "atacn_lead = get_dom_ASA_in_GM(snp_tab = LEDE_taba_f,encode_map = atac_tabn,encode_cids = c(7,8),cnames = c(\"ATAC_caQTL_ref\",\"ATAC_caQTL_alt\",\"ATAC_caQTL_more\",\"ATAC_caQTL_less\"))\n",
    "\n",
    "LEDE_ASAtab = get_dep_ASA_in_GM(snp_tab = atacn_lead,encode_map = atac_tabn,encode_cids = c(7,8),cnames = c(\"ATAC_dep_ref\",\"ATAC_dep_alt\",\"ATAC_dep_more\",\"ATAC_dep_less\"))\n",
    "\n",
    "print(head(LEDE_ASAtab))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2bc49e7",
   "metadata": {},
   "source": [
    "### add genotype annotation information in GM12878 for each SNP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "3025e088",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"***** Object of class 'vcfR' *****\"\n",
      "[1] \"***** Meta section *****\"\n",
      "[1] \"##fileformat=VCFv4.1\"\n",
      "[1] \"##FILTER=<ID=PASS,Description=\\\"All filters passed\\\">\"\n",
      "[1] \"##fileDate=20150218\"\n",
      "[1] \"##reference=ftp://ftp.1000genomes.ebi.ac.uk//vol1/ftp/technical/refe [Truncated]\"\n",
      "[1] \"##source=1000GenomesPhase3Pipeline\"\n",
      "[1] \"##contig=<ID=1,assembly=b37,length=249250621>\"\n",
      "[1] \"First 6 rows.\"\n",
      "[1] \n",
      "[1] \"***** Fixed section *****\"\n",
      "     CHROM POS       ID           REF ALT QUAL  FILTER\n",
      "[1,] \"1\"   \"565286\"  \"rs1578391\"  \"C\" \"T\" \"100\" \"PASS\"\n",
      "[2,] \"1\"   \"842362\"  \"rs28540380\" \"C\" \"T\" \"100\" \"PASS\"\n",
      "[3,] \"1\"   \"968226\"  \"rs3121560\"  \"C\" \"A\" \"100\" \"PASS\"\n",
      "[4,] \"1\"   \"1186502\" \"rs6603785\"  \"A\" \"T\" \"100\" \"PASS\"\n",
      "[5,] \"1\"   \"1353380\" \"rs11488462\" \"T\" \"C\" \"100\" \"PASS\"\n",
      "[6,] \"1\"   \"1355583\" \"rs1781133\"  \"C\" \"A\" \"100\" \"PASS\"\n",
      "[1] \n",
      "[1] \"***** Genotype section *****\"\n",
      "     FORMAT NA12878\n",
      "[1,] \"GT\"   \"0|1\"  \n",
      "[2,] \"GT\"   \"0|0\"  \n",
      "[3,] \"GT\"   \"0|0\"  \n",
      "[4,] \"GT\"   \"0|0\"  \n",
      "[5,] \"GT\"   \"0|0\"  \n",
      "[6,] \"GT\"   \"0|0\"  \n",
      "[1] \n",
      "[1] \"Unique GT formats:\"\n",
      "[1] \"GT\"\n",
      "[1] \n",
      "[1] \"\"\n",
      "  FORMAT NA12878\n",
      "1     GT     0|1\n",
      "2     GT     0|0\n",
      "3     GT     0|0\n",
      "4     GT     0|0\n",
      "5     GT     0|0\n",
      "6     GT     0|0\n",
      "  CHROM     POS         ID REF ALT QUAL FILTER\n",
      "1     1  565286  rs1578391   C   T  100   PASS\n",
      "2     1  842362 rs28540380   C   T  100   PASS\n",
      "3     1  968226  rs3121560   C   A  100   PASS\n",
      "4     1 1186502  rs6603785   A   T  100   PASS\n",
      "5     1 1353380 rs11488462   T   C  100   PASS\n",
      "6     1 1355583  rs1781133   C   A  100   PASS\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "vcf <- read.vcfR( \"/home/.../ALL.phase3_v5b.20130502.LeadSNPs.NA12878.vcf.gz\", verbose = FALSE )\n",
    "print(head(vcf))\n",
    "\n",
    "dfchrom = as.data.frame(vcf@fix)[,1:7]\n",
    "dfgt = as.data.frame(vcf@gt)[,1:2]\n",
    "\n",
    "print(head(dfgt))\n",
    "print(head(dfchrom))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "8194cc68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         Peak  Chr       Pos       RsID Ref Alt      AF Inside_Peak    P_Lead\n",
      "31476   29217 chr2  62408114 rs13004362   G   C 0.28005           1 0.2656618\n",
      "31476.1 29217 chr2  62408114 rs13004362   G   C 0.28005           1 0.2656618\n",
      "32257   29972 chr2  70237472 rs11126258   G   T 0.06500           1 0.3867709\n",
      "32257.1 29972 chr2  70237472 rs11126258   G   T 0.06500           1 0.3867709\n",
      "34401   32046 chr2 102372131  rs6543094   T   C 0.75000           1 0.3467704\n",
      "34401.1 32046 chr2 102372131  rs6543094   T   C 0.75000           1 0.3467704\n",
      "             Beta       SE  Log10_BF Chr  Pos_Left Pos_Right  Peak      GC\n",
      "31476    0.033513 0.046526 -1.720579   2  62407956  62408600 29217 43.7209\n",
      "31476.1  0.033513 0.046526 -1.720579   2  62407956  62408600 29217 43.7209\n",
      "32257   -0.023992 0.052718 -1.733415   2  70237165  70238157 29972 49.5468\n",
      "32257.1 -0.023992 0.052718 -1.733415   2  70237165  70238157 29972 49.5468\n",
      "34401    0.002010 0.028256 -2.047814   2 102371729 102372607 32046 43.1172\n",
      "34401.1  0.002010 0.028256 -2.047814   2 102371729 102372607 32046 43.1172\n",
      "        Peak_Height  P_caQTL          PMR P_eQTL_Coloc Peak_width Peak_center\n",
      "31476          2256 0.004783 0.0001908664 0.0017755589        644    62408278\n",
      "31476.1        2256 0.004783 0.0001908664 0.0017755589        644    62408278\n",
      "32257          1546 0.002939 0.0004331109 0.0011177076        992    70237661\n",
      "32257.1        1546 0.002939 0.0004331109 0.0011177076        992    70237661\n",
      "34401          3650 0.002201 0.0001069281 0.0002981116        878   102372168\n",
      "34401.1        3650 0.002201 0.0001069281 0.0002981116        878   102372168\n",
      "        Dist_SNP_to_peak_center SNP_centrality dep_Peak dep_Chr   dep_Pos\n",
      "31476                       164     0.25465839    29216    chr2  61880173\n",
      "31476.1                     164     0.25465839    29218    chr2  62409465\n",
      "32257                       189     0.19052419    29971    chr2  70229357\n",
      "32257.1                     189     0.19052419    29973    chr2  69744669\n",
      "34401                        37     0.04214123    32045    chr2 102608746\n",
      "34401.1                      37     0.04214123    32047    chr2 102374635\n",
      "           dep_RsID dep_Ref      dep_Alt  dep_AF dep_Inside_Peak dep_P_Lead\n",
      "31476    esv3591022       T        <CN0> 0.15495               0 0.16670031\n",
      "31476.1  rs71873378     TTG            T 0.72960               1 0.98765642\n",
      "32257     rs4852405       A            G 0.06500               1 0.65194383\n",
      "32257.1  esv3591185       T <INS:ME:ALU> 0.12755               0 0.04014632\n",
      "34401   rs201167739 GAGAGAC            G 0.05500               0 0.08103893\n",
      "34401.1  rs12712115       T            C 0.71000               1 0.94194529\n",
      "         dep_Beta   dep_SE dep_Log10_BF dep_Chr dep_Pos_Left dep_Pos_Right\n",
      "31476    0.082078 0.051806   -1.2503542       2     62377646      62378399\n",
      "31476.1  0.119447 0.044354   -0.3417915       2     62409069      62409622\n",
      "32257    0.060513 0.061095   -1.5032613       2     70229320      70230286\n",
      "32257.1  0.049525 0.038919   -1.5629444       2     70241969      70242785\n",
      "34401   -0.217807 0.053467    1.5393603       2    102370163     102370598\n",
      "34401.1  0.087004 0.033656   -0.5764270       2    102374371     102375039\n",
      "        dep_Peak  dep_GC dep_Peak_Height dep_P_caQTL      dep_PMR\n",
      "31476      29216 42.0424            1451    0.004138 3.388738e-05\n",
      "31476.1    29218 47.2924             417    0.003488 2.975831e-06\n",
      "32257      29971 45.7084             468    0.000511 3.659213e-08\n",
      "32257.1    29973 36.1077            2539    0.011073 1.390502e-03\n",
      "34401      32045 41.2844            3098    0.007009 1.326725e-03\n",
      "34401.1    32047 40.0598            3628    0.046021 4.255541e-03\n",
      "        dep_P_eQTL_Coloc dep_Peak_width dep_Peak_center\n",
      "31476       0.0008346143            753        62378022\n",
      "31476.1     0.0007336328            553        62409345\n",
      "32257       0.0002304469            966        70229803\n",
      "32257.1     0.0066066237            816        70242377\n",
      "34401       0.0009894685            435       102370381\n",
      "34401.1     0.0063812814            668       102374705\n",
      "        dep_Dist_SNP_to_peak_center dep_SNP_centrality LD_peak_dist\n",
      "31476                        497849        661.1540505        30256\n",
      "31476.1                         120          0.2169982         1067\n",
      "32257                           446          0.4616977         7858\n",
      "32257.1                      497708        609.9362745         4716\n",
      "34401                        238365        547.9655172         1787\n",
      "34401.1                          70          0.1047904         2537\n",
      "        enh_subcategory lead_anno_tss500 dep_anno_tss500 LEAD_TAD_id DEP_TAD_id\n",
      "31476         no-effect           Distal          Distal        <NA>       <NA>\n",
      "31476.1       no-effect           Distal          Distal        <NA>    tad_937\n",
      "32257         no-effect           Intron          Intron        <NA>       <NA>\n",
      "32257.1       no-effect           Intron          Intron        <NA>       <NA>\n",
      "34401         no-effect           Intron          Intron    tad_1000   tad_1000\n",
      "34401.1       no-effect           Intron          Intron    tad_1000   tad_1000\n",
      "        LEAD_AB DEP_AB caQTL_enh_GC_700bp dep_GC_700bp causal_label\n",
      "31476      <NA>   <NA>          0.4336662    0.4136947    no-effect\n",
      "31476.1    <NA>     A2          0.4336662    0.4607703    no-effect\n",
      "32257      <NA>   <NA>          0.5121255    0.4664765    no-effect\n",
      "32257.1    <NA>   <NA>          0.5121255    0.3680456    no-effect\n",
      "34401        A2     A2          0.4279601    0.3980029    no-effect\n",
      "34401.1      A2     A2          0.4279601    0.4065621    no-effect\n",
      "        enh_status_label LEAD_SNP_GMindex GM_RsID GM_REF GM_ALT GM_Pos\n",
      "31476          no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "31476.1        no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "32257          no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "32257.1        no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "34401          no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "34401.1        no-effect               NA    <NA>   <NA>   <NA>     NA\n",
      "        GM_phasing Dist_GM_SNP_to_peak_center lead_peak_id dep_peakid\n",
      "31476         <NA>                         NA        29217      29216\n",
      "31476.1       <NA>                         NA        29217      29218\n",
      "32257         <NA>                         NA        29972      29971\n",
      "32257.1       <NA>                         NA        29972      29973\n",
      "34401         <NA>                         NA        32046      32045\n",
      "34401.1       <NA>                         NA        32046      32047\n",
      "        depGMSNP_RsID depGMSNP_REF depGMSNP_ALT depGMSNP_Pos depGMSNP_phasing\n",
      "31476            <NA>         <NA>         <NA>           NA             <NA>\n",
      "31476.1          <NA>         <NA>         <NA>           NA             <NA>\n",
      "32257            <NA>         <NA>         <NA>           NA             <NA>\n",
      "32257.1          <NA>         <NA>         <NA>           NA             <NA>\n",
      "34401            <NA>         <NA>         <NA>           NA             <NA>\n",
      "34401.1          <NA>         <NA>         <NA>           NA             <NA>\n",
      "        depGMSNP_dist_to_dep_peak_center depGMSNP_refcount depGMSNP_altcount\n",
      "31476                                 NA                NA                NA\n",
      "31476.1                               NA                NA                NA\n",
      "32257                                 NA                NA                NA\n",
      "32257.1                               NA                NA                NA\n",
      "34401                                 NA                NA                NA\n",
      "34401.1                               NA                NA                NA\n",
      "        ATAC_caQTL_ref ATAC_caQTL_alt ATAC_caQTL_more ATAC_caQTL_less\n",
      "31476               NA             NA              NA              NA\n",
      "31476.1             NA             NA              NA              NA\n",
      "32257               NA             NA              NA              NA\n",
      "32257.1             NA             NA              NA              NA\n",
      "34401               NA             NA              NA              NA\n",
      "34401.1             NA             NA              NA              NA\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        ATAC_dep_ref ATAC_dep_alt ATAC_dep_more ATAC_dep_less GM_LEAD_RSID_geno\n",
      "31476             NA           NA            NA            NA               0|1\n",
      "31476.1           NA           NA            NA            NA               0|1\n",
      "32257             NA           NA            NA            NA               0|0\n",
      "32257.1           NA           NA            NA            NA               0|0\n",
      "34401             NA           NA            NA            NA               1|1\n",
      "34401.1           NA           NA            NA            NA               1|1\n",
      "        GM_LEAD_RSID_qual GM_LEAD_RSID_pass lead_geno_label lead_LCL_actstate\n",
      "31476                 100              PASS             HET               ALT\n",
      "31476.1               100              PASS             HET               ALT\n",
      "32257                 100              PASS             REF               REF\n",
      "32257.1               100              PASS             REF               REF\n",
      "34401                 100              PASS             ALT               ALT\n",
      "34401.1               100              PASS             ALT               ALT\n",
      "        lead_gm_actstate\n",
      "31476                het\n",
      "31476.1              het\n",
      "32257          homo_more\n",
      "32257.1        homo_more\n",
      "34401          homo_more\n",
      "34401.1        homo_more\n",
      "       \n",
      "           CM LOCAL no-effect\n",
      "  FALSE     0     0        12\n",
      "  TRUE   3908 16071     16128\n",
      "      \n",
      "          CM LOCAL no-effect\n",
      "  PASS  3908 16071     16128\n",
      "           \n",
      "              CM LOCAL no-effect\n",
      "  het        760  2735      2469\n",
      "  homo_less  631  2371      2636\n",
      "  homo_more  660  2519      2506\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "LEDE_ASAtab$GM_LEAD_RSID_geno = dfgt$NA12878[match(LEDE_ASAtab$RsID,dfchrom$ID)]\n",
    "LEDE_ASAtab$GM_LEAD_RSID_qual= dfchrom$QUAL[match(LEDE_ASAtab$RsID,dfchrom$ID)]\n",
    "LEDE_ASAtab$GM_LEAD_RSID_pass= dfchrom$FILTER[match(LEDE_ASAtab$RsID,dfchrom$ID)]\n",
    "\n",
    "LEDE_ASAtab$lead_geno_label = 'NA'\n",
    "LEDE_ASAtab$lead_geno_label[LEDE_ASAtab$GM_LEAD_RSID_geno=='0|1' | LEDE_ASAtab$GM_LEAD_RSID_geno=='1|0'] = 'HET'\n",
    "LEDE_ASAtab$lead_geno_label[LEDE_ASAtab$GM_LEAD_RSID_geno=='0|0' ] = 'REF'\n",
    "LEDE_ASAtab$lead_geno_label[LEDE_ASAtab$GM_LEAD_RSID_geno=='1|1' ] = 'ALT'\n",
    "\n",
    "LEDE_ASAtab$lead_LCL_actstate = c('REF','ALT')[c(1,2)[as.integer(LEDE_ASAtab$Beta>0)+1]]\n",
    "\n",
    "LEDE_ASAtab$lead_gm_actstate = 'NA'\n",
    "LEDE_ASAtab$lead_gm_actstate[LEDE_ASAtab$lead_geno_label=='HET'] = 'het'\n",
    "LEDE_ASAtab$lead_gm_actstate[(LEDE_ASAtab$lead_geno_label=='REF' & LEDE_ASAtab$lead_LCL_actstate=='REF') | (LEDE_ASAtab$lead_geno_label=='ALT' & LEDE_ASAtab$lead_LCL_actstate=='ALT')] = 'homo_more'\n",
    "LEDE_ASAtab$lead_gm_actstate[(LEDE_ASAtab$lead_geno_label=='REF' & LEDE_ASAtab$lead_LCL_actstate=='ALT') | (LEDE_ASAtab$lead_geno_label=='ALT' & LEDE_ASAtab$lead_LCL_actstate=='REF')] = 'homo_less'\n",
    "\n",
    "print(tail(LEDE_ASAtab))\n",
    "\n",
    "print(table(!is.na(LEDE_ASAtab$GM_LEAD_RSID_geno),LEDE_ASAtab$enh_subcategory))\n",
    "print(table(LEDE_ASAtab$GM_LEAD_RSID_pass,LEDE_ASAtab$enh_subcategory))\n",
    "print(table(LEDE_ASAtab$lead_gm_actstate[!duplicated(LEDE_ASAtab$Peak)],LEDE_ASAtab$enh_subcategory[!duplicated(LEDE_ASAtab$Peak)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c3b284f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "## write.csv(LEDE_ASAtab,'/home/.../Interim_Master_table_LEAD_DEP_pairs_noeffect_local_cm_AllspInfo.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python-3.9.5",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}