RNA_Seq_Data_Significant_Genes_msk8_paired.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FAIqfRcaRTt-"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1zkjQXQtgXRW"
   },
   "source": [
    "### Reads the data into a dataframe and renames the columns to remove unnecessary numbers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "E8JZsjDRRbla"
   },
   "outputs": [],
   "source": [
    "data_df = pd.read_excel(\"RNASeqData.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "L5SIC7ykeqLo"
   },
   "outputs": [],
   "source": [
    "data_df = data_df.rename(columns=lambda x: re.sub(\"w_[0-9]+_O\", \"w_O\",x))\n",
    "data_df = data_df.rename(columns=lambda x: re.sub(\"m_[0-9]+_O\", \"m_O\",x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "RfuwTmgERdju"
   },
   "outputs": [],
   "source": [
    "columns = list(data_df.columns) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "53N3U1QORlyV"
   },
   "outputs": [],
   "source": [
    "norm_cols = [x for x in columns if \"Norm\" in x] # get all normalized values\n",
    "healthy_control = [x for x in norm_cols if \"HC\" in x] # healthy control norm values\n",
    "norm_cols_base = [x for x in norm_cols if \"Base\" in x] # cf without treatment\n",
    "norm_cols_v2 = [x for x in norm_cols if \"V2\" in x] # cf with treatment\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qFE6PTePg0KZ"
   },
   "source": [
    "### Remove extra characters so we can keep patient identifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zzbolvnReQp1"
   },
   "outputs": [],
   "source": [
    "df_norm_base = data_df[norm_cols_base]\n",
    "df_norm_base = df_norm_base.rename(columns=lambda x: re.sub(\"Norm_Orkambi_0\", \"\",x))\n",
    "df_norm_base = df_norm_base.rename(columns=lambda x: re.sub(\"_Base\", \"\",x))\n",
    "df_norm_base = df_norm_base.sort_index(axis=1)\n",
    "df_norm_base = df_norm_base.drop(columns=['02'])\n",
    "\n",
    "df_norm_v2 = data_df[norm_cols_v2]\n",
    "df_norm_v2 = df_norm_v2.rename(columns=lambda x: re.sub(\"Norm_Orkambi_0\", \"\",x))\n",
    "df_norm_v2 = df_norm_v2.rename(columns=lambda x: re.sub(\"_V2\", \"\",x))\n",
    "df_norm_v2 = df_norm_v2.sort_index(axis=1)\n",
    "df_norm_v2 = df_norm_v2.drop(columns=['19'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nihjnYeuhCNC"
   },
   "source": [
    "### Function: get_mean_diff\n",
    "The function takes to lists and simply subtracts the corresponding values and finally takes the average of the differences.\n",
    "\n",
    "Making a function call can take extra time, so for simple operations like these, it may be helpful when used in conjunction with the **map()** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "AJ-A6XLhRoIp"
   },
   "outputs": [],
   "source": [
    "from statistics import mean\n",
    "def get_mean_diff(base, v2):\n",
    "  diff = [base[i] - v2[i] for i in range(len(base))] \n",
    "  mean_diff = mean(diff)\n",
    "  # mean of the vector \n",
    "  return mean_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "elvQrjV4mkfb"
   },
   "source": [
    "### Function: get_avg_diff_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hBTa-wkVRp0D"
   },
   "outputs": [],
   "source": [
    "import random\n",
    "from random import randrange\n",
    "from random import shuffle\n",
    "def get_avg_diff_count(expvalues):\n",
    "    num_shuffles=100000\n",
    "    base_copy = expvalues[0]\n",
    "    v2_copy = expvalues[1]\n",
    "    n_samples = len(base_copy)\n",
    "    #observed_avg_diff = get_mean_diff(base_copy, v2_copy)\n",
    "    obsdiff = [base_copy[i] - v2_copy[i] for i in range(n_samples)]\n",
    "    observed_avg_diff = mean(obsdiff)\n",
    "    count = 0\n",
    "    #avg_diffs = []\n",
    "    # try looping over shuffles\n",
    "    randomorder = np.array(random.choices([1,-1], k=n_samples*num_shuffles))\n",
    "    randomshape = (num_shuffles, n_samples)\n",
    "    randomorder = randomorder.reshape(randomshape)\n",
    "    randomdiff = np.multiply(randomorder, obsdiff)\n",
    "    avg_diffs = np.average(randomdiff, axis=1)\n",
    "\n",
    "    count=[]\n",
    "    count = np.count_nonzero(abs(avg_diffs) > abs(observed_avg_diff))\n",
    "    \n",
    "#    return count/num_shuffles, np.asarray(avg_diffs)\n",
    "    return count/num_shuffles\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "471OW3BMmw36"
   },
   "source": [
    "### Function: avg_diff_sig_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6cYRjRoLRumA"
   },
   "outputs": [],
   "source": [
    "def avg_diff_sig_test(df_norm_base, df_norm_v2):\n",
    "    fullzip = list(zip(df_norm_base.values, df_norm_v2.values))\n",
    "    all_p_sig_values = [get_avg_diff_count(fullzip[i]) for i in range(len(fullzip))]\n",
    "    #all_p_sig_values = map(get_avg_diff_count, fullzip)\n",
    "    return list(all_p_sig_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "Jj0uTrsFeQqr",
    "outputId": "6ae59e68-2b83-47ca-e466-d97eb731811e"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15570, 19)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_norm_v2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "PmpDSBStRwek",
    "outputId": "d9414b64-37ad-4434-a8d6-92e9a5f63372",
    "tags": [
     "outputPrepend"
    ]
   },
   "outputs": [],
   "source": [
    "norm_p_values = avg_diff_sig_test(df_norm_base, df_norm_v2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 164
    },
    "colab_type": "code",
    "id": "UjjjK_laeQq0",
    "outputId": "4d3ea549-fbb8-49de-f164-4b2921839883"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.28026, 0.28929, 0.0096, 0.6487, 0.10179]"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_p_values[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2quyhpBaRy0f"
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "pickle.dump( norm_p_values, open( \"norm_p_values.p\", \"wb\" ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "If3jPZW6R13a"
   },
   "source": [
    "**Bonferroni correction**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gb-bYne4R6TZ"
   },
   "source": [
    "In this method, we divide the threshold by the number of genes. So, if we take a default threshold of 5% (or 0.05), now our threshold would be `0.05/len(norm_p_values)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 170
    },
    "colab_type": "code",
    "id": "joWk5fWbSBRK",
    "outputId": "18fa802f-1bea-4efa-9c3b-873835adfc8e",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of significant genes:  0\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "norm_p_values = pickle.load( open( \"norm_p_values.p\", \"rb\" ) )\n",
    "count_significant_genes_bonferroni = 0\n",
    "for i, ele in enumerate(norm_p_values):\n",
    "  if ele < (0.05/len(norm_p_values)):\n",
    "    count_significant_genes_bonferroni += 1\n",
    "    print(i+1)\n",
    "print(\"Number of significant genes: \", count_significant_genes_bonferroni)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1wbHwRJkSH7K"
   },
   "source": [
    "**Benjamini–Hochberg procedure**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ltfby0LWSM4c"
   },
   "source": [
    " Put the individual P values in order, from smallest to largest. The smallest P value has a rank of i=1, then next smallest has i=2, etc. Compare each individual P value to its Benjamini-Hochberg critical value, (i/m)Q, where i is the rank, m is the total number of tests, and Q is the false discovery rate you choose. The largest P value that has P<(i/m)Q is significant, and all of the P values smaller than it are also significant, even the ones that aren't less than their Benjamini-Hochberg critical value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CSvKOcR5SIZJ"
   },
   "outputs": [],
   "source": [
    "def benjamin_hochberg_corrected_p_values(p_values, fdr):\n",
    "  num_genes = len(p_values)\n",
    "  import numpy as np\n",
    "  sorted_indices = np.argsort(p_values)\n",
    "  sorted_p_values = sorted(p_values)\n",
    "  benjamini_hochberg_critical_values = [((i+1)/num_genes)*fdr for i in range(num_genes)]\n",
    "  for i, p_value in reversed(list(enumerate(sorted_p_values))):\n",
    "    if p_value < benjamini_hochberg_critical_values[i]:\n",
    "      break\n",
    "  return i, sorted_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "fGZnls8PSRaf",
    "outputId": "9d878dd4-abbd-4a08-8268-5e79272e11a0",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of significant genes:  1\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "norm_p_values = pickle.load( open( \"norm_p_values.p\", \"rb\" ) )\n",
    "i, sorted_indices = benjamin_hochberg_corrected_p_values(norm_p_values, 0.25)\n",
    "print(\"Number of significant genes: \", i+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "post_drug.vs.pre_drug_MLELog2FC       -1.65495\n",
       "post_drug.vs.pre_drug_pVal         4.40299e-06\n",
       "post_drug.vs.pre_drug_pAdj           0.0356456\n",
       "control_Mean                             20.23\n",
       "pre-drug_Mean                            55.86\n",
       "post-drug_Mean                           17.08\n",
       "Name: 6444, dtype: object"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_df.loc[6444][14:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.08923"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_p_values[6444]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "01    26.220131\n",
       "04    38.436591\n",
       "06    10.362710\n",
       "07    11.041510\n",
       "09     5.162527\n",
       "10    10.750095\n",
       "12    18.467806\n",
       "13    11.509118\n",
       "14    18.359318\n",
       "15     7.410397\n",
       "16     5.962505\n",
       "17    10.882621\n",
       "18    67.777682\n",
       "20     3.313864\n",
       "21    12.973341\n",
       "22    12.691223\n",
       "24    32.851086\n",
       "25    13.435806\n",
       "26    11.839753\n",
       "Name: 6444, dtype: float64"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_norm_v2.loc[6444]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "01      6.437689\n",
       "04     35.964963\n",
       "06     64.415597\n",
       "07     14.532351\n",
       "09     19.714927\n",
       "10     13.925694\n",
       "12     34.883549\n",
       "13     24.114367\n",
       "14     27.847446\n",
       "15     29.548117\n",
       "16     94.710858\n",
       "17     18.286489\n",
       "18      5.678884\n",
       "20     17.214459\n",
       "21     10.864266\n",
       "22     15.517023\n",
       "24     47.494550\n",
       "25     20.203974\n",
       "26    578.546142\n",
       "Name: 6444, dtype: float64"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_norm_base.loc[6444]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "01     199.568365\n",
       "04      34.526364\n",
       "06      43.857428\n",
       "07      93.901346\n",
       "09      51.399631\n",
       "10      33.943880\n",
       "12     233.451446\n",
       "13      88.110189\n",
       "14      60.758065\n",
       "15      30.346715\n",
       "16      60.133878\n",
       "17     100.575692\n",
       "18     123.312905\n",
       "20      88.720672\n",
       "21      49.383029\n",
       "22    6757.663691\n",
       "24    3159.153584\n",
       "25      34.258913\n",
       "26     124.089284\n",
       "Name: 6659, dtype: float64"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_norm_base.loc[6659]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "01    541.088158\n",
       "04    183.312973\n",
       "06     71.387559\n",
       "07     41.011322\n",
       "09     79.650415\n",
       "10     19.019399\n",
       "12     25.854928\n",
       "13    247.969171\n",
       "14    134.635000\n",
       "15     59.283179\n",
       "16     68.142914\n",
       "17    133.700773\n",
       "18     28.953573\n",
       "20     85.055835\n",
       "21    140.710850\n",
       "22     21.649733\n",
       "24     55.664341\n",
       "25     49.609130\n",
       "26    159.298501\n",
       "Name: 6659, dtype: float64"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_norm_v2.loc[6659]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "nhaIX-MTSf2_",
    "outputId": "b092dc06-a8ec-484e-9ebe-e0c4d629affe"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[188,\n",
       " 1095,\n",
       " 10516,\n",
       " 13503,\n",
       " 5750,\n",
       " 11160,\n",
       " 10363,\n",
       " 361,\n",
       " 1549,\n",
       " 3290,\n",
       " 15332,\n",
       " 1392,\n",
       " 7204,\n",
       " 7868,\n",
       " 10593,\n",
       " 8384]"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "significant_genes = []\n",
    "for i in range(16):\n",
    "  significant_genes.append((sorted_indices[i]+1))\n",
    "\n",
    "\n",
    "significant_genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 170
    },
    "colab_type": "code",
    "id": "6wYK5OdSSirJ",
    "outputId": "37a1646e-69cb-421c-9898-466db0a45024",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       RowID            ID                                        Description\n",
      "187      188        ADAM22                    ADAM metallopeptidase domain 22\n",
      "360      361         AKIP1                     A-kinase interacting protein 1\n",
      "1094    1095          BCO2                          beta-carotene oxygenase 2\n",
      "1391    1392      C1orf115                chromosome 1 open reading frame 115\n",
      "1548    1549      CACNA2D2  calcium voltage-gated channel auxiliary subuni...\n",
      "3289    3290       EEF1DP3  eukaryotic translation elongation factor 1 del...\n",
      "5749    5750         KITLG                                         KIT ligand\n",
      "7203    7204  LOC105378698                                                NaN\n",
      "7867    7868        LRRC8D   leucine-rich repeat containing 8 family member D\n",
      "8383    8384       MIR3939                                      microRNA 3939\n",
      "10362  10363         PRDX1                                    peroxiredoxin 1\n",
      "10515  10516         PSMD1               proteasome 26S subunit, non-ATPase 1\n",
      "10592  10593        PTP4A2    protein tyrosine phosphatase type IVA, member 2\n",
      "11159  11160        RNF165                            ring finger protein 165\n",
      "13502  13503       TMEM181                          transmembrane protein 181\n",
      "15331  15332        ZNF593                            zinc finger protein 593\n"
     ]
    }
   ],
   "source": [
    "print(data_df.loc[data_df['RowID'].isin(significant_genes)][['RowID', 'ID', 'Description']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "jXcJvTCbTFNj",
    "outputId": "5bba2d0f-d2d2-42cc-f3d5-cc6b2b93788f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Norm_Orkambi_007_V2    2.201250e+07\n",
       "Norm_Orkambi_009_V2    2.599588e+07\n",
       "Norm_Orkambi_010_V2    2.801257e+07\n",
       "Norm_Orkambi_012_V2    2.576081e+07\n",
       "Norm_Orkambi_013_V2    2.615355e+07\n",
       "Norm_Orkambi_014_V2    2.786734e+07\n",
       "Norm_Orkambi_001_V2    3.455492e+07\n",
       "Norm_Orkambi_004_V2    3.376038e+07\n",
       "Norm_Orkambi_006_V2    2.666307e+07\n",
       "Norm_Orkambi_015_V2    2.896782e+07\n",
       "Norm_Orkambi_016_V2    2.332386e+07\n",
       "Norm_Orkambi_017_V2    2.449434e+07\n",
       "Norm_Orkambi_018_V2    2.312316e+07\n",
       "Norm_Orkambi_019_V2    2.542411e+07\n",
       "Norm_Orkambi_020_V2    2.495638e+07\n",
       "Norm_Orkambi_021_V2    2.707834e+07\n",
       "Norm_Orkambi_022_V2    2.569940e+07\n",
       "Norm_Orkambi_024_V2    2.229553e+07\n",
       "Norm_Orkambi_025_V2    3.108347e+07\n",
       "Norm_Orkambi_026_V2    2.474027e+07\n",
       "dtype: float64"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_df[norm_cols_v2].sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DCbRAr5yRp7A"
   },
   "outputs": [],
   "source": [
    "#draw histogram\n",
    "\n",
    "def draw_hist(d, observed):\n",
    "    print(observed)\n",
    "    hist,bin_edges = np.histogram(d,bins =100)    \n",
    "    plt.figure(figsize=[8,8])\n",
    "    plt.bar(bin_edges[:-1], hist, width = 0.01, color='#0504aa',alpha=0.7)\n",
    "    #plt.xlim(min(min(bin_edges),-4), max(max(bin_edges),4))\n",
    "    plt.grid(axis='y', alpha=0.75)\n",
    "    plt.xlabel('Value',fontsize=15)\n",
    "    plt.ylabel('Frequency',fontsize=15)\n",
    "    plt.xticks(fontsize=15)\n",
    "    plt.yticks(fontsize=15)\n",
    "    plt.ylabel('Frequency',fontsize=15)\n",
    "    plt.title('Normal Distribution Histogram',fontsize=15)\n",
    "    plt.axvline(x=observed, color='r', linestyle='dashed', linewidth=2)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "iFjKawc_tu-h"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAH8CAYAAABM96wGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA310lEQVR4nO3de7xcVX338c9XKkK4K0VSa4nGAlVa82i0moogaCwqklIV2/q0KhhtvLVSUClq8PIooRFbKNUUFG3rldKEqNwVCkWrROMNiZoKaEEUTILhIgrr+WPvYybDuc05M2dmzv68X695nczaa+/92/tMzv7NWmuvnVIKkiSpuR7U7wAkSVJ/mQxIktRwJgOSJDWcyYAkSQ1nMiBJUsOZDEiS1HAmAxoaSZYnKUkuHmXZeUmu6ENYHUtyaH0cB41TZ15dZ+R1Z5KNSf4tycGj1L8iyXkdxPCiJC/toP7yJLd1cgwdbHufevvz2sq7to9JxjHm/pK8tF62a/1+5PfzvElue8f6GBd0OWypK0wGNIwWJ3lSv4OYIX8DPBV4DvAO4GHAfyZ5W1u9ZcCbO9jui4CXdlD/bODZHdTvxD7A24B5beVfoTr2jT3a73TcQhXb1ZOsvyPVMS7oVUDSdPxavwOQOvRT4IfA3wJLur3xJDuXUu7u9nanYUMp5Yv1v68Ezk3ydmB5kitLKVcAlFKu68XOkzwYuL+U8kOq8z5jSil3AF+csGIflFJ+zoDG1moAP88aULYMaNgU4P8Bz0/yu+NVTLIgyeVJ7kqyqW5if3jL8pGm3j9L8pEkm4G1LeUvTvKhJHck+WGSl9TrnZjk5iQ/SXJqkge1bPPAJB9P8oN6v99K8letdbrgFOBm4FUt+92umyDJbyb5ZJIfJ7m77mJ4R73sXOCPgUNauiGWt24nydIkG4F7gN9o7yZo8RtJPl13Y9yU5FWtC0frvmhtjq+7Br5RL/r8SDzt9VrWnZPkH5L8KMk9Sb6cZPFo+0zyp0m+V//+Lkzym52c5PGM1k2Q5PlJ1tXnYlOS/05ySL34Z/XPD7Wc83n1ensn+XCS2+vPzBVJFrbt7yFJ/inJ5rreafXnqrTUGTlfz05yQZKtwJn1suPrc7Ulya1J1iZ5zBjn7WVJvp9ka5J/qff95CRfqsuuSPJb3TqXGgy2DGgYfYrqgvi3wItHq5Dk14ErgG8DfwrsCrwHuDTJwlLKvS3V/w44H3ghcF9L+anAv1FdOF8OfDjJ/wH2q98/EXgn8FXg4/U6jwA21Ov9jKpZ+BRgZ+DdUz/kbUop9yX5HPCAsQMtPlLvcymwGXg0cGC97B3AbwF7UnUvwPbf+v8AmA+8EbgL2DLOfs4B/gU4Azga+KckPyylfHqSh3ML8GdU5+vVVF0D4/ln4PnAScD3gFcAn0nyjFJKa5P97wO/ARxPdR7+HlhF1d0ykR2StP9tHDeZSzIfOK/ezwnATlSfj4fWVQ4DPkf1eflMXXZL/XM18BiqLqHb6vU/n+T/lFK+V9dZQdWtcxLVZ/pljPHZp/qdfAh4H1UyB/CbVInBjcDuVInkfyXZv5TS+vt9CrA38Fqqz8jpwN1U53MFcCfwD1Tn8g/HOycaMqUUX76G4gUsB26r//1Sqgv3/vX784ArWuq+h+oiuHtL2ZOpWhb+pH4/r37/H237GSn/UEvZ7sAvgO8CO7SUfwn4xBjxhirhPgn4n5byQ+vtHzTOsY7E8Lwxlr8buLvl/RXAeS3vtwJHjrP97c5X23buBvYd69y3HcOqtnqXAl8cK67Rjh84qH5/6AT1fge4H/iLljoPAr4JXNy2zy3AXi1lf1Vva+dxzsnI/sZ77Tra7wd4AXD7ONveta7/0rbyP6zLD2kp2wX4CfCB+v3D6t/JCW2frW8BZZT4T5/g/9EOVAnSz4A/bztvm4E9Wso+WW/z6S1ly+qyOVP9v+xr8F52E2hY/StwE2MPmnsycEmp+p0BKKV8CbgBeFpb3c8wustb1r2D6g/0laWU1taD71G1BgCQZKckpyT5HvBzqgTiXcCjRvm2OR2ZYPl64N2pRsF32qS7rpTyo0nW/Y+29+cDT0yyQ4f7nIwnUR33p0YKSin31+/bf6dfLqVsank/MqbiEUzsxfW+Wl+nTLDON4A96ub+xUl2mcR+oPqc/qSUcuVIQSnlTuDTbDum36VqabigpU4B1o6xzQd8npM8JcmlSW4HfknV4rMrsH9b1WvL9i0F3wPuZfuBkiOtFb8x/qFpmJgMaCiVUn5J1Wz5kiT7jVJlLnDrKOW3sq3ptrVsNJvb3t87RtlOLe9PpWruHWmSfhJV0zBt9abrEYwdN8AxwLVUzbw3Jlmf5PBJbnu87bb78Sjvf42qqbnb5gJbSyl3tZXfCsxJ8pCWss1tdUa6hSbzO/hWKeXa1hdVEjmmUsoG4Ciq7pjPArcl+WjdXTWeyXxO961//qStTvv71nV/pU4GL6FKpF5J1Q30JKrfVfv52Nz2/l7gZ3XS1VrGKOtqiJkMaJh9kOoP2htHWXYL1S1r7R5OdUdCq24+x/uFwBmllBWllMvqC8kvu7h96haGw4AvjFWnlPK/pZSXUjUxPxX4EXBBkodNYhednI/2c7wP1fGODDa8h+q2ulbtydhk3QLsmmROW/nDgbtKNcK/b0opnymlHEx1zo8Fnkk1lmI8k/mcjrTStCcWYyUa7b+/PwTmAEeVUs4rpVxD1XI01d+DZiGTAQ2t+o//31EN5pvbtvi/gWcn2W2kINXcBPOY/L3hU7EzVffAyD53YOyBXlP1Vqom2vdPVLGUcn+pbk08heqCMNKK0t6iMVV/NMr7dS1dKT9k28DFEc9qez/Zb5pfprrQvWCkIEnq9738nXaklLKllPJRqi6Ux9bFYx3jfwP7JHn6SEGd7DyXbcf0Daqk6qiWOgGOnGRIO1ONtWhNSl+EA8jVwg+Dht0HqAboLaK6D3/Ee4G/BC5Ocirb7ib4BvDvPYznUuDV9ZiBn1KNkH/I+KuM64BUt/TtCDyKKrH4Q2B5az9zqyR7ABdT3VHwnXr/x1N9w/x2Xe164KgkS6gu2DeXUm6eQnxHJHkX1bk/mupCf1TL8v8Ajk1yOlVf9jN44ORFN1ENkPuLJFuAX9QtKtsppXw7yceAM5Pszra7CQ6k+l33TZJXUrXAXER12+dvU7USfQSglHJvku8DL0ryTaqL+9dLKRcn+S/gE0neBNxO1c20M3Bave7tSf4ZOCXJL9h2N8HuTK4V53NUgwY/lOQc4HH1PjZ349g1O9gyoKFW9x+fPkr5T6guPPcAHwP+EbgKeFbZ/rbCbnttvZ9/pOrG+CbTu6Xw76i6Ay6mmsHup1Qju8cb0HYPVdLzeqpBZx+mGjC2uGybgOYsqn7kD1J94146xfiOA55AdXvc84BXl1JaB7p9hipZewFVYrAf1ch+WurcQ3VRfyJVUvHlcfb3ivp43gKsqbf3vLL9bYX98HWqZvv3Up3Xk6lug2ztwnoV1ViKy6iOcWQA3h9RJZHvoxoMGeCwsu22QoATgXOp7ur4GNW4gHOAO5hAKeUbVMnD71MNTPxTqkRlvFtG1TCpBqVKkoZJksuAB5dSDpmwsjQBuwkkacAleQbVN/uvAA+mulvkcKpv+NK0mQxI0uDbSvUsjjdTDUL8LtUERpN+UqU0HrsJJElqOAcQSpLUcCYDkiQ1XGPHDOy9995l3rx5/Q5jfD+pZxv99YlmNJUkaXzr1q27rZQy6gWlscnAvHnzuPbaB8xrIknSrJTkxrGW2U0gSVLDmQwMslWrqpckST3U2G6CofDKV1Y/l051plhJkiZmy4AkSQ1nMiBJUsOZDEiS1HAmA5IkNZzJgCRJDWcyIElSw3lr4SDziZKSpBlgy4AkSQ1nMiBJUsOZDAyyJz6xekmS1EOOGRhkX/lKvyOQJDWALQOSJDWcyYAkSQ1nMiBJUsOZDEiS1HAmA5IkNZx3E/TIkiVXb/d+9eqndb6RV7yiS9FIkjQ2k4FBtmpVvyOQJDXAjHYTJHlBkmuS3J7kniQbkpycZMeWOjckKW2vH42yrccmuTzJXUluTvL2JDvM5PFIkjQbzHTLwMOAzwOnAZuBJwPLgX2B17TU+yhwRsv7e1s3kmQv4DLgOuAoYD6wkiq5ObknkffDunXVT2chlCT10IwmA6WUD7QVfT7J7sCrk7y2lF89pu+WUsoXx9nUq4CdgaNLKXcAl9bbWZ5kRV02UKY0hmDhwuqnTy+UJPXQINxNcDuw44S1tncEcHHbRf/jVAnCId0KTJKkJuhLMpBkhyRzkjwNeB3wTy2tAgAvT3Jvki1JzkuyX9smDgSuby0opdwE3FUvkyRJk9SvuwnuBB5S//sjwAkty9YAXwR+CPwO8DbgqiS/W0rZUtfZi2rMQbtN9bIZ194NIEnSsOhXMrAImEM1gPCtwJnAMoBSyutb6l2V5BpgPfAy4H0ty0brSM8Y5dXCZCmwFGDu3LmsX79+qvE/wKJFWzuqP5l9L+igriRJU9WXZKCUMvJs3quT3AZ8OMnKUsrGUep+M8kG4AktxZuAPUfZ9B6M3mIwsq1VwCqAhQsXlgULFkwp/tEsX95Zy8CJJ05+392MU5KkdoMwgHAkMXjUBPVav/FfT9vYgCSPBHahbSyBJEka3yAkA39Q//z+aAuTHAQcAKxrKb4QeHaS3VrKjgHuBq7sRZB9ce211UuSpB6a0W6CJBdRTRb0LeA+qkTgeOATpZSNSZ4LvAT4NHAz1bf/k4GbgHNbNvV+qrsQzk9yKvBoqsmL3juIcwxMmZMNSZJmwEyPGfgy8FJgHvBL4H+AN1Nd3AF+AOxDNVBwT6o5CC4CTmq9yJdSNiU5nGrg4VqqcQKnUyUEkiSpAzM9A+FbgLeMs/zrwOGT3NZ1wGFdCm0wLV1a/fSBRZKkHhqEMQMayz//c/WSJKmHTAYkSWo4kwFJkhrOZECSpIYzGZAkqeFMBiRJarh+PahIk/GEJ0xcR5KkaTIZGGTr1k1cR5KkabKbQJKkhjMZkCSp4UwGBllSvSRJ6iGTAUmSGs4BhANiyZKrt3u/evXT+hSJJKlpbBmQJKnhTAYkSWo4kwFJkhrOZECSpIZzAOEg+8AH+h2BJKkBTAYG2dKl/Y5AktQAdhNIktRwJgODbNWq6iVJUg/ZTTDIXvnK6qfdBZKkHrJlQJKkhjMZkCSp4UwGJElqOJMBSZIazmRAkqSGMxmQJKnhvLVwkJXS7wgkSQ1gy4AkSQ1nMiBJUsOZDAyyJz6xekmS1EOOGRhkX/lKvyOQJDWALQOSJDWcyYAkSQ1nMiBJUsOZDEiS1HAmA5IkNZx3EwyyV7yi3xFIkhrAZKBPliy5euJKq1b1PhBJUuPZTSBJUsOZDAyydeuqlyRJPWQ3wSBbuLD66dMLJUk9ZDIwJNrHGKxe/bQ+RSJJmm3sJpAkqeFmNBlI8oIk1yS5Pck9STYkOTnJji11kuSkJD9IcneS/0yyYJRtPTbJ5UnuSnJzkrcn2WEmj0eSpNlgplsGHgZ8HjgOOAL4IPC3wHtb6rwJeAtwKnAksBW4LMm+IxWS7AVcBhTgKODtwPHAKb0/hJnR2i0wqdsQJUmaohkdM1BK+UBb0eeT7A68OslrgYdQJQPvLqWcCZDkC8ANwGuAk+v1XgXsDBxdSrkDuLTezvIkK+oySZI0CYMwZuB2YKSbYBGwO/DJkYWllDuBtVQtCSOOAC5uu+h/nCpBOKSn0UqSNMv05W6Cum//IcATgNcB/1RKKUkOBO4Dvtu2yreBY1reHwh8rrVCKeWmJHfVy9b2KvaZdPwhZ/c7BElSA/Tr1sI7qZIBgI8AJ9T/3gvYWkq5r63+JmBOkh1LKffW9TaPst1N9bJZYeOeB/Q7BElSA/QrGVgEzAGeDLwVOBNYVi8bbYadjLJsrHpjztCTZCmwFGDu3LmsX7++o6DHs2jR1q5tazK6Gbskqdn6kgyUUr5S//PqJLcBH06ykuqb/W5JdmhrHdgTuKuU8ov6/aa6rN0ejN5iMLLfVcAqgIULF5YFCxZM4yi2t3x590f8L1u/AoCzFpz4gGUnnrig6/uTJDXTIAwgHEkMHgVcD+wAPKatzoH1shHX12W/kuSRwC5t9Yba4hvXsvjGWTH8QZI0wAYhGfiD+uf3gWuAO4AXjixMModqvoELW9a5EHh2kt1ayo4B7gau7Gm0kiTNMjPaTZDkIqrJgr5FddfAH1BNFvSJUsrGus57gLck2UT1Lf8NVEnLGS2bej/VXQjnJzkVeDSwHHivcwxIktSZmR4z8GXgpcA84JfA/wBvprq4j3gP1cX/zVQzFl4LPKuUcutIhVLKpiSHUw08XEs1TuB0qoRAkiR1YKZnIHwL1VTD49UpwLvq13j1rgMO6150kiQ10yCMGZAkSX3Ur3kGNAkb99i/3yFIkhrAZGCAHX/oOf0OQZLUAHYTSJLUcCYDkiQ1nMnAAFu95mBWrzm432FIkmY5kwFJkhrOAYSzxJIlD3xQ0urVT+tDJJKkYWPLgCRJDWcyIElSw5kMSJLUcCYDkiQ1nAMIB9hZjz+h3yFIkhrAZGCAXTLv+f0OQZLUAHYTSJLUcCYDA2zxDRew+IYL+h2GJGmWs5tggC372mmA3QWSpN6yZUCSpIYzGZAkqeFMBiRJajiTAUmSGs5kQJKkhjMZkCSp4by1cIAtOeqqfocgSWoAk4FZbMmSq7d7v3r10/oUiSRpkNlNIElSw5kMDLCVVxzLyiuO7XcYkqRZzm6CATZ/y3f6HYIkqQFMBoZU+3gASZKmym4CSZIazmRAkqSGMxmQJKnhTAYkSWo4BxAOsEv2O7LfIUiSGsBkYICdteDEfocgSWoAuwkkSWo4k4EBNn/zBuZv3tDvMCRJs5zdBANs5ZXHAT69UJLUW7YMSJLUcCYDkiQ1nMmAJEkNZzIgSVLDmQxIktRwJgOSJDWctxYOsOMPObvfIUiSGmBGWwaSvDDJBUn+N8nWJOuS/ElbnRuSlLbXj0bZ1mOTXJ7kriQ3J3l7kh1m7mh6b+OeB7BxzwP6HYYkaZab6ZaBNwDfB/4auA14DvDRJHuXUs5oqfdRoPX9va0bSbIXcBlwHXAUMB9YSZXcnNyz6CVJmoVmOhk4spRyW8v7zyX5DaokofXif0sp5YvjbOdVwM7A0aWUO4BLk+wOLE+yoi4besvWrwB698CiJUuu3u796tVP68l+JEmDbUa7CdoSgRFfBfbpcFNHABe3XfQ/TpUgHDLF8AbO4hvXsvjGtf0OQ5I0yw3C3QSLqJr7W708yb1JtiQ5L8l+bcsPBK5vLSil3ATcVS+TJEmT1Ne7CZIcTtXn//KW4jXAF4EfAr8DvA24KsnvllK21HX2AjaPsslN9TJNgd0GktRMfUsGksyjGii4ppRy7kh5KeX1LdWuSnINsB54GfC+lmVltM2OUT6yz6XAUoC5c+eyfv36KcU+mkWLtnZtW7+yprvbbj/eibbbzfMjSRpcfUkGkjwUuBC4CXjJeHVLKd9MsgF4QkvxJmDPUarvwegtBiPbWgWsAli4cGFZsGBBJ2GPa/nyqyeuNEXXXLNrV7Zz4okLtns/Uczt9SVJs9OMjxlIMgf4NLAj8NxSyp2TXLX1G//1tI0NSPJIYBfaxhJIkqTxzfSkQ78GfAr4beCIUsqPJ7HOQcABwLqW4guBZyfZraXsGOBu4MruRdxfG/fYn4177N/vMCRJs9xMdxOcRTXR0OuBhyZ5SsuyrwLPpOo2+DRwM9W3/5OpuhPOban7fuB1wPlJTgUeDSwH3jtb5hgAOP7Qc/odgiSpAWY6GVhc//z7UZY9CvgB1ZwD76MaE3A7cBFwUutFvpSyqb4T4UxgLdU4gdOpEgJJktSBGU0GSinzJlHt8Elu6zrgsGkFJEmSBmLSIY1h9ZqDWb3m4H6HIUma5XyEcYO0TyokSRLYMiBJUuOZDEiS1HAmA5IkNZzJgCRJDWcyIElSw3k3wQA76/En9DsESVIDmAwMsEvmPb/fIUiSGsBuAkmSGs6WgQG2+IYLgMFpIWiftGj16qf1KRJJUjeZDAywZV87DRicZECSNDvZTSBJUsOZDEiS1HAmA5IkNZzJgCRJDWcyIElSw5kMSJLUcN5aOMCWHHVVv0OQJDWALQOSJDWcyYAkSQ1nMjDAVl5xLCuvOLbfYUiSZjnHDAyw+Vu+0+8QJEkNYMuAJEkNZzIgSVLDmQxIktRwJgOSJDWcyYAkSQ3n3QQD7JL9jux3CJKkBjAZGGBnLTix3yGMa8mSq7d7v3r10/oUiSRpOjrqJkhyUK8CkSRJ/dHpmIGvJ/lykr9MsmcvAtI28zdvYP7mDf0OQ5I0y3WaDBwOXAesAG5O8rEkz0qS7oemlVcex8orj+t3GJKkWa6jZKCU8vlSyl8A+wKvAR4BXAzcmOQdSeb3IEZJktRDU7q1sJRyZynlg6WUpwP7AzcAJwHfSXJlkj/qYoySJKmHpjzPQJJ5SZYDlwBPBT4LLAVuBT6R5PSuRChJknqqo1sLk8wBXgC8DDgY+D7wz8C5pZRb6mrnJHkZ8PfAX3cxVg246d5q6K2KktQfnc4z8CNgB+B84JmllCvGqPdl4PZpxKUB0H5xliTNTp0mA28C/q2UsmW8SqWUbwKPmnJUkiRpxnSUDJRSzupVIHqg4w85u98hSJIaoNMxAx8EdimlHDPKso8Bd5ZSvDG+SzbueUC/Q5AkNUCndxM8CzhvjGX/DiyeXjiSJGmmdZoM/Drw0zGWbQL2mV44arVs/QqWrV/R7zAkSbNcp8nAjcDTx1j2dOCH0wtHrRbfuJbFN67tdxiSpFmu02TgXOCNSV6dZFeAJLsmWQacCDjiTZKkIdNpMnAq8K/AGcCWJHcAW4AzgQ/Xy8eU5IVJLkjyv0m2JlmX5E/a6iTJSUl+kOTuJP+ZZMEo23psksuT3JXk5iRvT7JDh8cjSVLjdXpr4f3AcUlOAw4DHko1udDnSinfmcQm3kA1a+FfA7cBzwE+mmTvUsoZdZ03AW8BTgCur9e5LMlBpZQfASTZC7iM6gmKRwHzgZVUyc3JnRyTJElN1+mkQwCUUjYAG6aw6pGllNta3n8uyW9QXfDPSLITVTLw7lLKmQBJvkD1IKTXsO1C/ypgZ+DoUsodwKVJdgeWJ1lRl0mSpEmYUjKQZH/gN4Gd2peVUj471npticCIr1J9uwdYBOwOfLJlnTuTrAWOYFsycARwcdtF/+NU3RSHAI66mwV8VoEkzYxOJx16LPAJ4LFARqlSqJ5d0IlFVM39AAcC9wHfbavzbaB1oqMDgc9tt+NSbkpyV71sViQDG/fYv98hSJIaoNOWgQ8AOwJHU13A753OzpMcTtUq8PK6aC9gaynlvraqm4A5SXYspdxb19s8yiY31ctmheMPPaffIUiSGqDTZOD/AC8upXx6ujtOMg/4KLCmlHJuy6IyWvVRlo1Vb7TykX0uBZYCzJ07l/Xr108+4AksWrS1a9uaLVasuGi794sX7ztu/YnOYTd/X5KkbTpNBjYyyjiBTiV5KHAhcBPwkpZFm4DdkuzQ1jqwJ3BXKeUXLfX2HGXTezB6iwEApZRVwCqAhQsXlgULFkztAEaxfLmP+53IiScuGHf5ROdwovUlSVPT6TwDxwMnJXn0VHeYZA7waaruhueWUu5sWXw91ZiDx7StdmC9rLXegW3bfSSwS1u9obZ6zcGsXnNwv8OQJM1ynbYMvBt4BHB9khsY5Vt4KeXJY62c5NeATwG/DfxBKeXHbVWuAe4AXgi8s15nDnAk9Tf62oXACUl2K6X8rC47BrgbuLLDY5IkqdE6TQa+Wb+m6iyqiYZeDzw0yVNaln21lHJPkvcAb0myiW2TDj2IatbDEe8HXgecn+RU4NHAcuC9zjGgXvJ2R0mzUaczEL5smvsbecTx34+y7FFUkwu9h+ri/2bgYcC1wLNKKbe2xLGpvhPhTKrbCDcDp1MlBJIkqQNTnXQoVJMOPRL4Wlu//5hKKfMmUacA76pf49W7jmpKZA0Jv1VL0mDqdAAh9RMK/5fqccZXAQfU5ecn+auuRidJknqu0xkITwDeQTXt7+fZfhbAK4A/Ad7XpdikaWlviQBbIyRpNJ12E7waeGspZcUojwveADh/bhed9fgT+h2CJKkBOk0G9gXWjbHsfrowIZG2uWTe8/sdgqbJcRKShkGnycD3qJ4KePkoy57OtgcOSV030YV1tG4BSdLEOk0G3gecleRe4Ly6bJ8kx1LNB/CKLsbWeItvuACwhUCS1FudzjNwdpK9gLcCp9TFnwXuApaXUj7a5fgabdnXTgNMBiRJvdXxPAOllNOSvB9YRDUp0E+BL5RStnQ7OM1uNutL0mCY0qRD9fMALu5yLJIkqQ86nWdg2UR1SilnTT0cSZI00zptGThznGWl/mkyIEnSEOl0AOEDpi9OsifwbOCNVDMQSo3hPAKSZoMpjRloVUrZDHwiyR7AB4BDp7tNSZI0c6adDLT4PrCwi9trvCVHXdXvECRJDdDxUwtHk2QucDxVQiBJkoZIp3cT/IRtAwVH7AjsBtwDHN2luCRJ0gzptJvgH3lgMnAP8EPgolLK7V2JSgCsvOJYAI4/9Jw+RzJ7OOBPkh6o07sJlvcoDo1i/pbv9DsESVIDdGXMgCRJGl6djhn4Pg/sJhhTKeXRHUckSZJmVKdjBs4DXgzMAS4FfgzsAzwLuBP4RFejkyRJPddpMrAJ2Ag8t5Ry50hhkl2BTwNbSinv7GJ80oxygKGkJup0zMCrgdNaEwGAUspW4O/q5ZIkaYh02jKwB/DwMZbtC+w6vXDU6pL9jux3CJKkBug0GbgAOC3JHcDaUsrPkzwEeD5wKrC22wE22VkLTux3CJKkBug0GfhL4Fzgk0BJ8jOq2QdDlSj8ZVejkyRJPdfppENbgD9K8jjgSVRdBj8CvlxKua4H8TXa/M0bANi45wF9jkSSNJtN6amFpZRvAd/qcixqs/LK4wCfXihJ6q2OZyBMsk+SU5NcnmRD3UpAktcneWr3Q5QkSb3UUTKQ5MnAd4E/Bm4AHgM8pF488hhjSZI0RDrtJjgd+DzVo4ofBLysZdmXgD/tUlzSjGifZEiSmqjTZOAJwFGllPuTpG3Z7VRTE0szwgu5JHVHp2MGtgC/PsayRwO3Ti8cSZI00zpNBtYApyRpfRphSbI38DfA+V2LTJIkzYhOuwneBFwOXAesq8veTzWQ8PvAW7sXmo4/5Ox+hyBJaoBOJx3alOQpwP8FDqd6bPFPgbOBj5RSft79EJvLyYYkSTNh0slAkp2ophz+f6WUc4BzehaVpK7wkcySJmPSyUAp5Z4kTwJ26GE8arFs/QrABxbNJt2+OHuxl9QNnQ4gvABY0oM4NIrFN65l8Y0+CFKS1FudDiC8mOoRxnOBz1LdSlhaK5RSPtul2KS+cy4DSU3QaTLwr/XPo+tXu4LdCFLfmLxImooJk4EklwCvLaVsAB4FhOpOgv8Gftbb8CS18mIvqRcm0zLwTGAPgFLKjUl2AFYBTyql3NjL4CRJUu91/AjjWvtzCSRJ0pDqdMyAZtDGPfbvdwiSpAaYbDJQJlk2oSSPAU4AngIcBFxVSjm0rc4NwH5tq95aStm3rd5jgTOApwKbqWZCPKWUct9UYhs0xx/qvE7qLuclkDSaySYDFyf5ZVvZ5aOUUUqZ6DHGjwOeA3wR2HGceh+lutCPuLd1YZK9gMuonpNwFDAfWEnV9XHyBDFIkqTaZJKBU7q8z7WllDUASc4D9h6j3i2llC+Os51XATsDR5dS7gAuTbI7sDzJirpMkiRNYMJkoJTS1WSglHJ/lzZ1BHBx20X/48CpwCHA0E/dt3rNwQAsOeqqPkciSZrNpno3wUx4eZJ7k2xJcl6S9jEEBwLXtxaUUm4C7qqXSZKkSRjUuwnWUI0p+CHwO8DbgKuS/G4pZUtdZy+qQYPtNtXLJEnSJAxkMlBKeX3L26uSXAOsB14GvK+16iirZ4xykiwFlgLMnTuX9evXdyHayqJFW7u2rV9Z08NtqyfaP1MT/e46/Qx2+7PQzf8DkobXQCYD7Uop30yyAXhCS/EmYM9Rqu/B6C0GlFJWUc2eyMKFC8uCBQu6FuPy5b2bJvaaa3bt2bbVXSeeuGC79xN9LtrrT6Tbn7NO9y9pdhrkMQOjaf3Gfz1tYwOSPBLYhbaxBJIkaWxDkQwkOQg4AFjXUnwh8Owku7WUHQPcDVw5g+FJkjTUZrybIMkcqkmHAB4B7J7kBfX7zwLPAF4CfBq4merb/8nATcC5LZt6P/A64PwkpwKPBpYD750tcwyc9fgT+h2CJKkB+jFmYB/gU21lI+8fBfygrvM+qjEBtwMXASe1XuRLKZuSHA6cSTWnwGbgdKqEYFa4ZN7z+x2CJKkBZjwZKKXcwMRPPTx8ktu6DjhsujFJktRkQ3E3QVMtvuECwBaCJmt/sJAk9YLJwABb9rXTAJOB2cynCEoaBCYDkrZjgiI1j8mA1EU260saRiYDkvrKlgip/4Zi0iFJktQ7JgOSJDWcyYAkSQ3nmIEBtuSoq/odgmaYAxAl9YMtA5IkNZwtA5I0y3nHhiZiMjDAVl5xLADHH3pOnyPRbDUT3RJeiKTBZzIwwOZv+U6/Q5AkNYBjBiRJajiTAUmSGs5uAkmN4hgG6YFMBiSNy4unNPvZTSBJUsPZMjDALtnvyH6HIE3IlgNp+JkMDLCzFpzY7xAkSQ1gN4EkSQ1nMjDA5m/ewPzNG/odhiRplrObYICtvPI4wKcXSpJ6y2RA0kBzgKLUeyYDktQhExTNNiYDkrpqJp6EKKm7TAYkDTW/pUvT590EkiQ1nMmAJEkNZzfBADv+kLP7HYLUdb0eU2C3QfeN9jvzvM4uJgMDbOOeB/Q7BElSA9hNIElSw9kyMMCWrV8B+MAiSeOza0TTZTIwwBbfuBYwGZD6zbkTps+EZbDZTSBJUsPZMiCpI8P2LXnY4pX6wZYBSZIazmRAkqSGs5tAkvqs24Pr+tE14gDB4WYyMMA27rF/v0OQGmcQxxh4oVWvmQwMsOMPPaffIUgDZxAv1t02DFM29/v30OkxTFR/usuHnWMGJElqOFsGJDVav7/havYaptYEk4EBtnrNwQAsOeqqPkcizRwvzmqqfiYPJgOSpL4bpm/Rs9GMJwNJHgOcADwFOAi4qpRyaFudAG8G/hLYG/gy8LpSyvq2eo8FzgCeCmwGzgZOKaXc19ODkCTNKFuMeqsfLQOPA54DfBHYcYw6bwLeQpU0XA+8AbgsyUGllB8BJNkLuAy4DjgKmA+spBoUeXIvD0CSNFhsWZieftxNsLaU8shSyguBb7UvTLITVTLw7lLKmaWUy4AXAgV4TUvVVwE7A0eXUi4tpbwfOAV4Q5Lde34UkiTNEjPeMlBKuX+CKouA3YFPtqxzZ5K1wBFs+9Z/BHBxKeWOlnU/DpwKHAKs7VrQkjQOv5U+0HSb9Wdjt8Agf04GcZ6BA4H7gO+2lX+7XtZa7/rWCqWUm4C72upJkqRxDOLdBHsBW0cZBLgJmJNkx1LKvXW9zaOsv6leNvTOevwJ/Q5BknpiNn7zH2aDmAxANT6gXUZZNla90cpJshRYCjB37lzWr18/jRC3t2jR1q5ta8TWRYdV26b725bUOytWXLTd+8WL9x23fqd/P9r/dk33789Efwt78fet1yY6R91ePpqJzlv752TRou2Xd/MaNZGUMup1c2Z2npwH7N16a2GSZcA/AA9pbR1IcgKwvJSyS/3+x8A/llJOadvmVqrbC08bb98LFy4s1157bdeOxSxX0lgmmud+0LbXbhj/vnX72QNT6e/v9e+lU0nWlVIWjrZsEFsGrgd2AB4DbGgpbx8jcD1tYwOSPBLYpa3e0Fp8wwUAXDLv+X2ORNJsNsgD2/plGBOg6RjEAYTXAHdQ3U4IQJI5wJHAhS31LgSenWS3lrJjgLuBK2cgzp5b9rXTWPa1cRs4JEmatn7MQDiHatIhgEcAuyd5Qf3+s6WUu5K8B3hLkk1sm3ToQVSzDY54P/A64PwkpwKPBpYD72273VCS1IGmfStWf7oJ9gE+1VY28v5RwA3Ae6gu/m8GHgZcCzyrlHLryAqllE1JDgfOpJpTYDNwOlVCIEmSJqkfkw7dwLY7A8aqU4B31a/x6l0HHNa14CSpB/ymrUE3iGMGJEnSDBrEuwkkSZoWW2M6Y8uAJEkNZ8vAAFty1FX9DkGS1AAmA5IkdWi0bohhnqzJbgJJkhrOZGCArbziWFZecWy/w5AkzXJ2Ewyw+Vu+0+8QJEkNYMuAJEkNZzIgSVLDmQxIktRwJgOSJDWcyYAkSQ3n3QQD7JL9jux3CJKkBjAZGGBnLTix3yFIkhrAZECShoxP5BtMw/x7cczAAJu/eQPzN2/odxiSpFnOloEBtvLK4wCfXihJ6i1bBiRJajiTAUmSGs5kQJKkhnPMgCSpcYZ55H8v2DIgSVLDmQxIktRwdhMMsOMPObvfIUiSGsBkYIBt3POAfocgSWoAuwkkSWo4k4EBtmz9CpatX9HvMCRJs5zJwABbfONaFt+4tt9hSJJmOZMBSZIazmRAkqSGMxmQJKnhTAYkSWo4kwFJkhrOSYcG2MY99u93CJKkBjAZGGDHH3pOv0OQJDWA3QSSJDWcyYAkSQ1nMjDAVq85mNVrDu53GJKkWc5kQJKkhjMZkCSp4UwGJElqOJMBSZIazmRAkqSGMxmQJKnhBjIZSPLSJGWU16ta6iTJSUl+kOTuJP+ZZEEfw+66sx5/Amc9/oR+hyFJmuUGfTriw4C7W97/T8u/3wS8BTgBuB54A3BZkoNKKT+auRB755J5z+93CJKkBhj0ZODLpZSt7YVJdqJKBt5dSjmzLvsCcAPwGuDkmQxSkqRhNpDdBJOwCNgd+ORIQSnlTmAtcES/guq2xTdcwOIbLuh3GJKkWW7Qk4GNSX6ZZEOSV7aUHwjcB3y3rf6362WzwrKvncayr53W7zAkSbPcoHYT3EI1HuBLwA7AnwDvTzKnlHI6sBewtZRyX9t6m4A5SXYspdw7oxFLkjSkBjIZKKVcDFzcUnRhkocAJyf5+5Fqo6yacZaRZCmwFGDu3LmsX7++OwEDixY9YGjD9K3p4bYlSQOtm9eoiQxkMjCG84AXAfOoWgB2S7JDW+vAnsBdpZRfjLaBUsoqYBXAwoULy4IFC7oW3PLlV3dtW+2uuWbXnm1bkjSYTjxxwYzta9DHDIymUN1KuAPwmLZlB9bLJEnSJA1TMvDHwG3AjcA1wB3AC0cWJpkDHAlc2JfoJEkaUgPZTZDk36kGD36dqgXgmPr1ulLK/cA9Sd4DvCXJJrZNOvQg4Iz+RC1J0nAayGQA2AC8HHgk1aDA64A/L6X8S0ud91Bd/N8MPAy4FnhWKeXWGY61Z5YcdVW/Q5AkNcBAJgOllJOAkyaoU4B31S9JkjRFwzRmQJIk9YDJwABbecWxrLzi2H6HIUma5Qaym0CV+Vu+0+8QJEkNYMuAJEkNZzIgSVLDmQxIktRwJgOSJDWcyYAkSQ3n3QQD7JL9jux3CJKkBjAZGGBnLTix3yFIkhrAbgJJkhrOZGCAzd+8gfmbN/Q7DEnSLGc3wQBbeeVxgE8vlCT1li0DkiQ1nMmAJEkNZzIgSVLDmQxIktRwJgOSJDWcyYAkSQ3nrYUD7PhDzu53CJKkBjAZGGAb9zyg3yFIkhrAbgJJkhrOZGCALVu/gmXrV/Q7DEnSLGcyMMAW37iWxTeu7XcYkqRZzmRAkqSGMxmQJKnhTAYkSWo4kwFJkhrOZECSpIZz0qEBtnGP/fsdgiSpAUwGBtjxh57T7xAkSQ1gN4EkSQ1nMiBJUsOZDAyw1WsOZvWag/sdhiRpljMZkCSp4UwGJElqOJMBSZIazmRAkqSGMxmQJKnhTAYkSWo4ZyAcYGc9/oR+hyBJagCTgQF2ybzn9zsESVID2E0gSVLDmQwMsMU3XMDiGy7odxiSpFnOboIBtuxrpwF2F0iSemuoWwaSPDbJ5UnuSnJzkrcn2aHfcUmSNEyGtmUgyV7AZcB1wFHAfGAlVYJzch9DkyRpqAxtMgC8CtgZOLqUcgdwaZLdgeVJVtRlkiRpAsPcTXAEcHHbRf/jVAnCIf0JSZKk4TPMycCBwPWtBaWUm4C76mWSJGkShjkZ2AvYPEr5pnqZJEmahGEeMwBQRinLGOUkWQosrd9uTbKhBzHtDdzWjQ1l5B9rDu7G5oZJ185hg3kOu8PzOH2ewynKry4CXTuH+421YJiTgU3AnqOU78HoLQaUUlYBq3oXEiS5tpSysJf7mO08h9PnOewOz+P0eQ6nbybO4TB3E1xP29iAJI8EdqFtLIEkSRrbMCcDFwLPTrJbS9kxwN3Alf0JSZKk4TPMycD7gZ8D5yd5Zj0eYDnw3j7PMdDTboiG8BxOn+ewOzyP0+c5nL6en8OUMupYu6GQ5LHAmcBTqcYJnA0sL6Xc18+4JEkaJkOdDEiSpOkb5m6CGTXVhyIl2SPJh5JsSrIlyb8ledhMxDxopnIOkzypPn/fq9fbkORtSXaaqbgHyXQfzpXkQUnWJSlJntfLWAfZdM5jkqOTfDnJ3UluT3JRkl16HfOgmcbfxIVJLqnP3U+TXJbk92ci5kGT5DFJPpDka0nuS3LFJNfr+nVlmG8tnDHTfCjSJ4ADgOOA+4FTgdVAoyYPmMY5PKaueyrwXeD3gHfUP/+4hyEPnC49nOs44BE9CXBITOc8JjmOqmtyBXAC1QRnh9Gwv6VTPYf1HV+XAV8B/rwuPgG4JMnvlVJu7GXcA+hxwHOALwI7drBe968rpRRfE7yAN1PNa7B7S9mJVFMf7z7Oek+lmgDp6S1lT67Lntnv4xqSc/jro5Qtrc/hfv0+rmE4hy119wJ+Ahxbn7/n9fuYhuk8Uk388jPgFf0+hn6/pnEOXwXcB+zZUrZXXfaX/T6uPpzHB7X8+zzgikms05Prit0EkzPVhyIdAdxaSvnPkYJSypeA79fLmmRK57CU8pNRir9a/9yne+ENhek+nOsdwH8Bl/cgtmEy1fP4ovrnh3sV2BCZ6jl8MPBLYGtL2da6LKOuMYuVUu6fwmo9ua6YDEzOVB+K9ID1at+eYL3ZqJsPllpE1TTWi+mkB9mUz2GS3wNeBvxNz6IbHlM9j79P9Zk7NskPk/wiyX8nWdS7UAfWVM/hv9d1VibZJ8k+wOlUrQyf6lGss01PrismA5Mz1Yci+TClbbpyLpLsC/wt8C+lv/NJ9MN0zuEZwD+WUr7X7aCG0FTP475U/bQnA28EjgTuBC5K8vAuxzjopnQOSyk3A8+gGu9za/06Gnj2GK2AeqCeXFdMBiavo4cidWG92Wha5yLJjsAnqZoV/7qLcQ2Tjs9hkhdTXcTe2aughtBUPosPAnYFji2l/Fsp5SJgCVV/92u6HuHgm8pncS5V3/g6qibtI+p/fybJb/UiyFmq69cVk4HJ6fihSBOst+cE681GUz2HACQJ8BHq0bellE3dDG5IdHwOkzwYOI1qtPGDkuwJ7F4v3qVtOu+mmOpn8af1zytGCurWqXXAY7sT2tCY6jk8gerOixeUUi6qE6o/pkqo7MKanJ5cV0wGJmeqD0V6wHq1sfp8ZrPpPljqdKpbmI4qpTTt3I2YyjncBfhN4L1Uf0Q2AV+rl32cbYMxm2Sqn8VvU33zah/oFqoxLE0y1XN4IPCtUsovRgpKKfcC36K6PVET68l1xWRgcqb6UKQLgX2TPG2kIMlC4NH1siaZ8oOlkrwZeC3wklLK1b0LceBN5RxupeqjbX39Sb3sJODPehPqQJvqZ/HTVBf+Z4wUJNkDeCLbEqymmOo5vBE4qO7yAyDJQ4CDgBt6EOds1JvrSr/vsxyGF9WgjFuAS4FnUt3nvhV4Z1u97wHntJVdBPwP1SCZJVSjka/q9zENyzkE/pTq29iHgKe0vR4wB8Fsfk3nc9i2fB7NnmdgOv+fV9fr/gXwXKoL30+Avfp9XMNwDqkSp18An6nP3/OoLmC/AB7f7+Pqw3mcA7ygfn2BqoVk5P2c0c5hXdb160rfT8awvKj6BD9HlfneQnXP9g5tdW4Azm0r27O+kG0G7gA+Cuzd7+MZlnMInFtfuEZ7vbTfxzQM53CUbTQ6GZjOeaQaQPhPwO31upcBv9vv4xmyc3g48J9UYzB+SpVQHdrv4+nTORz5vzjaa94457Dr1xUfVCRJUsM5ZkCSpIYzGZAkqeFMBiRJajiTAUmSGs5kQJKkhjMZkCSp4UwGJI0pyaeTfGOc5Wcm2VTPIjfedg5NUpIc1P0oJU2XyYCk8XyMavrYx7UvSLID1Uxp55dSfj7jkUnqGpMBSeNZA9wFvHiUZc8AHk6VMEgaYiYDksZUStlK9YCeY0ZZ/GLgVuCWJB9P8oMkdyX5VpK/SjLm35ck8+pug+e1lZ+b5Nq2soOSfCbJz+rXp5Ls24XDk1QzGZA0kY8Bv53kiSMFSR4M/BHwSWBfqgelLAOeA/wzcArwxunuOMljgP8CdgL+L/BS4HHA2iTtjxKWNEW/1u8AJA28C6keiPJiYF1d9mzgocDHSilfAC4HqC/QV1M9je0VwLunue+3AT8CjijVc+9J8nWq57Y/h+rpd5KmyZYBSeOqBwf+B/Cilm/jx1A9m/6LSXZKckqS7wE/p3oc7buARyWZ7heOZ9b7vj/Jr9Xb+z7Vk9wWTnPbkmomA5Im42PAbwFPTbITcBRVq0ABTgX+BlhF9W39ScA76/V2muZ+96bqbvhF2+vRwCOnuW1JNbsJJE3G56gGC74YmAvsxra7CF4InFFKWTFSOclzJ9jePfXPHdvKH9r2/qdULQNnj7KN2yYOW9JkmAxImlAp5b4kn6K68D8C+HYp5ev14p2pugeAX80/MNqtiK1+TPUN/3da1tsVeCpV98OIy4GDgHV1K4SkHjAZkDRZHwNeQ3UXwVtbyi8FXl2PGfgp8Gpg3BkJSyn3J1kD/HWSG6kGKB4P3N1WdTnwJeAzST5I1RrwCOBZwLmllCumd0iSwDEDkibvC1QD9wJ8vKX8tcBVwD8CHwS+yeTuIngN1W2DZ9XrfoyqO+JXSinfAZ5CNfHRKqo7G06haon43pSPRNJ2YsubJEnNZsuAJEkNZzIgSVLDmQxIktRwJgOSJDWcyYAkSQ1nMiBJUsOZDEiS1HAmA5IkNZzJgCRJDff/Ac7zjV/YJDHAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_hist(norm_p_values, 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_sorted_indices = np.argsort(norm_p_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  187,  1094, 10515, 13502,  5749])"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_sorted_indices[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_sorted_p_values = sorted(norm_p_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1e-05,\n",
       " 0.00017,\n",
       " 0.00022,\n",
       " 0.00026,\n",
       " 0.00026,\n",
       " 0.00029,\n",
       " 0.00029,\n",
       " 0.00031,\n",
       " 0.00033,\n",
       " 0.00038]"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_sorted_p_values[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_critical_values = [((i+1)/15000)*0.1 for i in range(15000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[6.6666666666666675e-06,\n",
       " 1.3333333333333335e-05,\n",
       " 2e-05,\n",
       " 2.666666666666667e-05,\n",
       " 3.3333333333333335e-05,\n",
       " 4e-05,\n",
       " 4.666666666666667e-05,\n",
       " 5.333333333333334e-05,\n",
       " 5.9999999999999995e-05,\n",
       " 6.666666666666667e-05]"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_critical_values[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Akternative method\n",
    "story \n",
    "https://puolival.github.io/multipy/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from multipy.fwer import sidak\n",
    "\n",
    "#sidak_pvals = sidak(norm_p_values, alpha=.10)\n",
    "#print(sidak_pvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from multipy.fdr import qvalue\n",
    "#_, qvals = qvalue(norm_p_values, threshold=0.10)\n",
    "#len(qvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sum(qvals<0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "#tmp_sorted_qvals = sorted(qvals)\n",
    "#tmp_sorted_qvals[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.stats.multitest as ssm\n",
    "norm_p_values_np = np.array(norm_p_values, float)\n",
    "\n",
    "rej, s_adjpval = ssm.multipletests(norm_p_values_np, alpha=0.05, method=\"s\")[:2]\n",
    "rej, fdr_adjpval = ssm.multipletests(norm_p_values_np, alpha=0.05, method=\"fdr_i\")[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.14418476962410787,\n",
       " 0.9291453732841632,\n",
       " 0.967475995077026,\n",
       " 0.9825554213421549,\n",
       " 0.9825554213421549,\n",
       " 0.9890668393337582,\n",
       " 0.9890668393337582,\n",
       " 0.9919930734499568,\n",
       " 0.9941361445966126,\n",
       " 0.9973086825283916]"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_sorted_s_adjpval = sorted(s_adjpval)\n",
    "tmp_sorted_s_adjpval[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.1557,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507,\n",
       " 0.4235855895196507]"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_sorted_fdr_adjpval = sorted(fdr_adjpval)\n",
    "tmp_sorted_fdr_adjpval[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000010</td>\n",
       "      <td>0.000170</td>\n",
       "      <td>0.000220</td>\n",
       "      <td>0.000260</td>\n",
       "      <td>0.000260</td>\n",
       "      <td>0.000290</td>\n",
       "      <td>0.000290</td>\n",
       "      <td>0.000310</td>\n",
       "      <td>0.000330</td>\n",
       "      <td>0.000380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.155700</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "      <td>0.423586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.144185</td>\n",
       "      <td>0.929145</td>\n",
       "      <td>0.967476</td>\n",
       "      <td>0.982555</td>\n",
       "      <td>0.982555</td>\n",
       "      <td>0.989067</td>\n",
       "      <td>0.989067</td>\n",
       "      <td>0.991993</td>\n",
       "      <td>0.994136</td>\n",
       "      <td>0.997309</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0         1         2         3         4         5         6  \\\n",
       "0  0.000010  0.000170  0.000220  0.000260  0.000260  0.000290  0.000290   \n",
       "1  0.155700  0.423586  0.423586  0.423586  0.423586  0.423586  0.423586   \n",
       "2  0.144185  0.929145  0.967476  0.982555  0.982555  0.989067  0.989067   \n",
       "\n",
       "          7         8         9  \n",
       "0  0.000310  0.000330  0.000380  \n",
       "1  0.423586  0.423586  0.423586  \n",
       "2  0.991993  0.994136  0.997309  "
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame([tmp_sorted_p_values[:10],tmp_sorted_fdr_adjpval[:10],tmp_sorted_s_adjpval[:10]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "another_qvalue = estimate(norm_p_values_np, pi0=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.085635,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168,\n",
       " 0.10628096793349168]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_sorted_another_qvalue = sorted(another_qvalue)\n",
    "tmp_sorted_another_qvalue[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(np.array(tmp_sorted_another_qvalue)<0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# slightly modified code from \n",
    "# https://github.com/nfusi/qvalue/blob/master/qvalue/qvalue.py\n",
    "# changes appear to be transition from python2 to python3\n",
    "\n",
    "import scipy as sp\n",
    "def estimate(pv, m=None, verbose=False, lowmem=False, pi0=None):\n",
    "    \"\"\"\n",
    "    Estimates q-values from p-values\n",
    "\n",
    "    Args\n",
    "    =====\n",
    "\n",
    "    m: number of tests. If not specified m = pv.size\n",
    "    verbose: print verbose messages? (default False)\n",
    "    lowmem: use memory-efficient in-place algorithm\n",
    "    pi0: if None, it's estimated as suggested in Storey and Tibshirani, 2003.\n",
    "         For most GWAS this is not necessary, since pi0 is extremely likely to be\n",
    "         1\n",
    "\n",
    "    \"\"\"\n",
    "    assert(pv.min() >= 0 and pv.max() <= 1), \"p-values should be between 0 and 1\"\n",
    "\n",
    "    original_shape = pv.shape\n",
    "    pv = pv.ravel()  # flattens the array in place, more efficient than flatten()\n",
    "\n",
    "    if m is None:\n",
    "        m = float(len(pv))\n",
    "    else:\n",
    "        # the user has supplied an m\n",
    "        m *= 1.0\n",
    "\n",
    "    # if the number of hypotheses is small, just set pi0 to 1\n",
    "    if len(pv) < 100 and pi0 is None:\n",
    "        pi0 = 1.0\n",
    "    elif pi0 is not None:\n",
    "        pi0 = pi0\n",
    "    else:\n",
    "        # evaluate pi0 for different lambdas\n",
    "        pi0 = []\n",
    "        lam = np.arange(0, 0.90, 0.01)\n",
    "        counts = np.array([(pv > i).sum() for i in np.arange(0, 0.9, 0.01)])\n",
    "        for l in range(len(lam)):\n",
    "            pi0.append(counts[l]/(m*(1-lam[l])))\n",
    "\n",
    "        pi0 = np.array(pi0)\n",
    "\n",
    "        # fit natural cubic spline\n",
    "        tck = sp.interpolate.splrep(lam, pi0, k=3)\n",
    "        pi0 = sp.interpolate.splev(lam[-1], tck)\n",
    "        if verbose:\n",
    "            print(\"qvalues pi0=%.3f, estimated proportion of null features \" % pi0)\n",
    "\n",
    "        if pi0 > 1:\n",
    "            if verbose:\n",
    "                print(\"got pi0 > 1 (%.3f) while estimating qvalues, setting it to 1\" % pi0)\n",
    "            pi0 = 1.0\n",
    "\n",
    "    assert(pi0 >= 0 and pi0 <= 1), \"pi0 is not between 0 and 1: %f\" % pi0\n",
    "\n",
    "    if lowmem:\n",
    "        # low memory version, only uses 1 pv and 1 qv matrices\n",
    "        qv = np.zeros((len(pv),))\n",
    "        last_pv = pv.argmax()\n",
    "        qv[last_pv] = (pi0*pv[last_pv]*m)/float(m)\n",
    "        pv[last_pv] = -np.inf\n",
    "        prev_qv = last_pv\n",
    "        for i in range(int(len(pv))-2, -1, -1):\n",
    "            cur_max = pv.argmax()\n",
    "            qv_i = (pi0*m*pv[cur_max]/float(i+1))\n",
    "            pv[cur_max] = -np.inf\n",
    "            qv_i1 = prev_qv\n",
    "            qv[cur_max] = min(qv_i, qv_i1)\n",
    "            prev_qv = qv[cur_max]\n",
    "\n",
    "    else:\n",
    "        p_ordered = np.argsort(pv)\n",
    "        pv = pv[p_ordered]\n",
    "        qv = pi0 * m/len(pv) * pv\n",
    "        qv[-1] = min(qv[-1], 1.0)\n",
    "\n",
    "        for i in range(len(pv)-2, -1, -1):\n",
    "            qv[i] = min(pi0*m*pv[i]/(i+1.0), qv[i+1])\n",
    "\n",
    "        # reorder qvalues\n",
    "        qv_temp = qv.copy()\n",
    "        qv = np.zeros_like(qv)\n",
    "        qv[p_ordered] = qv_temp\n",
    "\n",
    "    # reshape qvalues\n",
    "    qv = qv.reshape(original_shape)\n",
    "\n",
    "    return qv\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15570"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(norm_p_values_np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., ..., 0., 0., 0.])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(tmp_sorted_another_qvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.500277"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(norm_p_values_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Healthy control vs CF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "RNA_Seq_Data_Significant_Genes.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}