6_remove_drug_interference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Remove-durg-interference-(Remove-Drug-comfoundings)\" data-toc-modified-id=\"Remove-durg-interference-(Remove-Drug-comfoundings)-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Remove durg interference (Remove Drug comfoundings)</a></span></li><li><span><a href=\"#Load-data\" data-toc-modified-id=\"Load-data-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Load data</a></span></li><li><span><a href=\"#Drug-AE-pairs-significantly-associated-with-pandemic\" data-toc-modified-id=\"Drug-AE-pairs-significantly-associated-with-pandemic-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Drug-AE pairs significantly associated with pandemic</a></span></li><li><span><a href=\"#Drug-significantly-associated-with-AE-during-pandemic\" data-toc-modified-id=\"Drug-significantly-associated-with-AE-during-pandemic-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Drug significantly associated with AE during pandemic</a></span></li><li><span><a href=\"#AEs-that-satisfy-both-significance-check\" data-toc-modified-id=\"AEs-that-satisfy-both-significance-check-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>AEs that satisfy both significance check</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove durg interference (Remove Drug comfoundings)\n",
    "\n",
    "\n",
    "For each top-N SE, we find the corresponding drugs, then count the occurance and non-occurance in 2019 and 2020. After that, calculate ROR, p-value, then use Fisher's test tojudge sig/insig, and multiple test to correct the  p-value.\n",
    "\n",
    "We check two aspects.First, the adverse event (such as hallucination) should be significantly associated with the therapy of at least one drug (like Pimavanserin). Second, the formed drug-adverse event pair (like Pimavanserin-hallucination) should be significantly associated with the pandemic\n",
    "\n",
    "**Exchange the order of step 3 and step 4 will not change the results.**\n",
    "\n",
    "**This section is computational expensive as containing nest loops.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "from tqdm import tqdm\n",
    "from collections import Counter\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pickle\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "# %matplotlib notebook\n",
    "pd.set_option('display.max_columns', None)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from sklearn.preprocessing import Normalizer, StandardScaler, MinMaxScaler\n",
    "import faiss\n",
    "\n",
    "# calculate 95% confidential interval\n",
    "def weird_division(n, d):\n",
    "    return n / d if d else 0\n",
    "\n",
    "def CI(ROR, A, B, C, D):\n",
    "    ror = np.log(ROR)\n",
    "    sq = 1.96*np.sqrt(weird_division(1, A) + weird_division(1, B) + weird_division(1, C) +weird_division(1, D))\n",
    "    CI_up = np.exp(ror + sq)\n",
    "    CI_down = np.exp(ror - sq)\n",
    "    return CI_up, CI_down\n",
    "def format_tex(float_number):\n",
    "    exponent = np.floor(np.log10(float_number))\n",
    "    mantissa = float_number/10**exponent\n",
    "    mantissa_format = str(mantissa)[0:3]\n",
    "    if float_number!=0:\n",
    "        return \"$< {0}\\times10^{{{1}}}$\".format(mantissa_format, str(int(exponent)))\n",
    "    else:\n",
    "        return \"$< 0 \\times10^{0}$\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_pd_US_pro = pickle.load(open('../Data/pandemic/all_pd_US_pro.pk', 'rb'))\n",
    "# drug_code_dic = pickle.load(open('../Data/parsed/drug_code_dic.pk', 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this MeDRA_dic, key is string of PT_name, value is a list:\n",
    "# [PT, PT_name, HLT,HLT_name,HLGT,HLGT_name,SOC,SOC_name,SOC_abbr]\n",
    "# medra_se_disease_dic = pickle.load(open('../Data/curated/AE_dic.pk', 'rb'))\n",
    "MedDRA_dic_all = pickle.load(open('../Data/curated/AE_mapping.pk', 'rb'))\n",
    "drug_dic = pickle.load(open('../Data/curated/drug_mapping.pk', 'rb'))\n",
    "# SE_dic = pickle.load(open('../Data/curated/AE_dic.pk', 'rb'))\n",
    "\n",
    "drug_dic_pd = pd.DataFrame(drug_dic, index=['Drugbank_ID', 'code'])\n",
    "drug_dic_pd = drug_dic_pd.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Drug-AE pairs significantly associated with pandemic\n",
    "Some SE of the Top-AE list are missed in the final dataframe, it's because this AE doesn't have significant drug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load condition SE_uncondition_2019_sig_over pop uncondition\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 150/150 [00:13<00:00, 11.26it/s]\n",
      "100%|██████████| 150/150 [00:15<00:00,  9.45it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_uncondition_2019_sig_over\n",
      "load condition SE_uncondition_2019_sig_under pop uncondition\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 33/33 [00:03<00:00,  9.63it/s]\n",
      "100%|██████████| 33/33 [00:03<00:00, 10.84it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_uncondition_2019_sig_under\n",
      "load condition SE_male_2019_sig_over pop male\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 42/42 [00:01<00:00, 38.69it/s]\n",
      "100%|██████████| 42/42 [00:01<00:00, 26.65it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_male_2019_sig_over\n",
      "load condition SE_male_2019_sig_under pop male\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 7/7 [00:00<00:00, 16.13it/s]\n",
      "100%|██████████| 7/7 [00:00<00:00, 17.24it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_male_2019_sig_under\n",
      "load condition SE_female_2019_sig_over pop female\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 87/87 [00:06<00:00, 13.04it/s]\n",
      "100%|██████████| 87/87 [00:08<00:00, 10.12it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_female_2019_sig_over\n",
      "load condition SE_female_2019_sig_under pop female\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 17/17 [00:01<00:00, 10.45it/s]\n",
      "100%|██████████| 17/17 [00:01<00:00, 11.36it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_female_2019_sig_under\n",
      "load condition SE_young_2019_sig_over pop young\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 18/18 [00:00<00:00, 134.77it/s]\n",
      "100%|██████████| 18/18 [00:00<00:00, 90.55it/s] \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_young_2019_sig_over\n",
      "load condition SE_young_2019_sig_under pop young\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 3/3 [00:00<00:00, 77.06it/s]\n",
      "100%|██████████| 3/3 [00:00<00:00, 178.79it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_young_2019_sig_under\n",
      "load condition SE_adult_2019_sig_over pop adult\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 96/96 [00:04<00:00, 20.87it/s]\n",
      "100%|██████████| 96/96 [00:05<00:00, 16.67it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_adult_2019_sig_over\n",
      "load condition SE_adult_2019_sig_under pop adult\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 8/8 [00:00<00:00, 44.62it/s]\n",
      "100%|██████████| 8/8 [00:00<00:00, 55.96it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_adult_2019_sig_under\n",
      "load condition SE_elderly_2019_sig_over pop elderly\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 52/52 [00:01<00:00, 42.24it/s]\n",
      "100%|██████████| 52/52 [00:01<00:00, 30.51it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_elderly_2019_sig_over\n",
      "load condition SE_elderly_2019_sig_under pop elderly\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 6/6 [00:00<00:00, 12.91it/s]\n",
      "100%|██████████| 6/6 [00:00<00:00, 14.04it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved SE_elderly_2019_sig_under\n"
     ]
    }
   ],
   "source": [
    "## Contingency table:\n",
    "# A: 2020_A, the # of drug-SE pair in 2020\n",
    "# B: 2020_B , the # of not this certain drug-SE pair in 2020\n",
    "# C:2019_A\n",
    "# D: 2019_B\n",
    "\n",
    "condition_list = ['SE_uncondition_2019_sig_over'\n",
    ", 'SE_uncondition_2019_sig_under', 'SE_male_2019_sig_over', 'SE_male_2019_sig_under',\n",
    "                 'SE_female_2019_sig_over', 'SE_female_2019_sig_under', \n",
    "                 'SE_young_2019_sig_over', 'SE_young_2019_sig_under', 'SE_adult_2019_sig_over', 'SE_adult_2019_sig_under',\n",
    "                 'SE_elderly_2019_sig_over', 'SE_elderly_2019_sig_under']\n",
    "\n",
    "for condition in condition_list:\n",
    "    \n",
    "    locals()[condition] = pickle.load(open('../Data/pandemic/results/'+condition+'_step2.pk', 'rb'))   \n",
    "    pop = condition.split('_')[1]\n",
    "    \n",
    "    print('load condition', condition, 'pop', pop)\n",
    "    df = locals()[condition]\n",
    "    if len(df) ==0:\n",
    "        print('there is no significant SE in this situation. ')\n",
    "        continue\n",
    "    \n",
    "    df = df.sort_values('p_corrected', ascending=True)\n",
    "    top_SE_list = list(df.SE)  # find the all the SE\n",
    "\n",
    "\n",
    "    # for top_SE in top_SE_list:\n",
    "    SE_list = top_SE_list\n",
    "    yr_list = [2019, 2020]\n",
    "\n",
    "    \"\"\"Find the corresponding Drugs of each Top SE\"\"\"\n",
    "    ind = ['2019-03-11'<i<'2019-09-31' or '2020-03-11'<i<'2020-09-31' for i in all_pd_US_pro['receipt_date']]\n",
    "    reports_19_20= all_pd_US_pro[ind]\n",
    "\n",
    "    top_drug_dic = {}\n",
    "    for i in range(len(SE_list)):\n",
    "        se = SE_list[i]\n",
    "        indx_SE = [se in j for j in reports_19_20.SE]\n",
    "        drug_list = reports_19_20[indx_SE].drugs  # find all the drugs which co-occured with the SE\n",
    "        drug_list = list(itertools.chain(*list(drug_list)))\n",
    "        drug_list = list(set(drug_list))  # remove duplicated drugs\n",
    "        top_drug_dic[i] = drug_list\n",
    "\n",
    "    \"\"\"initialize the Data frame, count the occurance for each SE-drug pairs \"\"\"\n",
    "    se_matrix = pd.DataFrame()  #{'SE': SE_list, 'name':list(SE_dic_df.index), 'medra_ID': list(SE_dic_df['medra_ID'])}\n",
    "\n",
    "    for yr in yr_list: # for a certain year\n",
    "        st = str(yr)+'-03-10'\n",
    "        end = str(yr)+'-09-31'\n",
    "        ind = [st<i<end for i in all_pd_US_pro['receipt_date']]\n",
    "        input_df= all_pd_US_pro[ind]\n",
    "        n_report = len(input_df)\n",
    "        ### split age_grop\n",
    "        input_df['age'] = pd.to_numeric(input_df['age'], errors='coerce')\n",
    "        bins = [1, 20, 65, max(input_df.age)+1]\n",
    "        age_labels = ['Young', 'Adult','Elderly']\n",
    "        input_df['age_group'] = pd.cut(input_df.age, bins, right = False, labels= age_labels)\n",
    "\n",
    "\n",
    "        \"\"\"Limit the reports into specific population!!!!!!!!!!!\"\"\"\n",
    "        if pop =='male':\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df[input_df.gender=='1']\n",
    "        elif pop =='female':\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df[input_df.gender=='2']\n",
    "        elif pop =='young':\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df[input_df.age_group=='Young']\n",
    "        elif pop =='adult':\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df[input_df.age_group=='Adult']\n",
    "        elif pop =='elderly':\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df[input_df.age_group=='Elderly']\n",
    "        else:\n",
    "            locals()['all_pd_US_pro_'+str(yr)] = input_df\n",
    "\n",
    "\n",
    "        se_, drug_, rank = [], [], []\n",
    "        A =[]\n",
    "        for i in tqdm(range(len(SE_list))): # dive into each SE\n",
    "            se = SE_list[i]    \n",
    "            name = locals()['all_pd_US_pro_'+str(yr)]\n",
    "            indx_SE = [se in j for j in name.SE]\n",
    "            drug_list = top_drug_dic[i]\n",
    "\n",
    "            #  dive into each SE-drug pair\n",
    "            SE_reports = name[indx_SE]\n",
    "            for drug in drug_list: # check the \n",
    "                index_SE_drug = [drug in j for j in SE_reports.drugs]        \n",
    "\n",
    "                se_.append(se) # record the se\n",
    "                drug_.append(drug) # record the drug\n",
    "                n_A = sum(index_SE_drug)      \n",
    "                A.append(n_A)  # record A\n",
    "                rank.append(i+1)\n",
    "        B = [n_report - i for i in A]\n",
    "        se_matrix['SE_rank'] = rank\n",
    "        se_matrix['SE'] = se_\n",
    "        se_matrix['drug'] = drug_\n",
    "        se_matrix[str(yr)+'_A'] = A\n",
    "        se_matrix[str(yr)+'_B'] = B\n",
    "        \n",
    "        \n",
    "\n",
    "    # insert the the drug name and SE name to the dataframe    \n",
    "    se_matrix['SE_name'] = se_matrix.apply(lambda row: list(MedDRA_dic_all[MedDRA_dic_all.PT==row.SE]['PT_name'])[0], axis=1)\n",
    "    se_matrix['drug_name'] = se_matrix.apply(lambda row: list(drug_dic_pd[drug_dic_pd['Drugbank_ID']==row.drug].index)[0], axis=1)\n",
    "\n",
    "    \n",
    "    se_matrix['2019_Delta'] = (se_matrix['2020_A'] - se_matrix['2019_A'])/se_matrix['2019_A']\n",
    "    # se_matrix\n",
    "\n",
    "    \"\"\"Calculate the Fisher's test, and correct the p-value\"\"\"\n",
    "    # Fisher's test\n",
    "    se_matrix_2019 = se_matrix\n",
    "    se_matrix_2019['2019_ROR'] = se_matrix_2019.apply(lambda row: stats.fisher_exact([[row['2020_A'], row['2020_B']], [row['2019_A'], row['2019_B']]])[0], axis = 1)\n",
    "    se_matrix_2019['2019_ROR_CI_upper'] = se_matrix_2019.apply(lambda row: CI(row['2019_ROR'],row['2020_A'], row['2020_B'],  row['2019_A'], row['2019_B'])[0], axis = 1)\n",
    "    se_matrix_2019['2019_ROR_CI_lower'] = se_matrix_2019.apply(lambda row: CI(row['2019_ROR'], row['2020_A'], row['2020_B'], row['2019_A'], row['2019_B'])[1], axis = 1)\n",
    "    \n",
    "    se_matrix_2019['p_value'] = se_matrix_2019.apply(lambda row: stats.fisher_exact([[row['2020_A'], row['2020_B']], [row['2019_A'], row['2019_B']]])[1], axis = 1)\n",
    "\n",
    "    # multipletests\n",
    "    se_matrix_2019['sig'], se_matrix_2019['p_corrected']  = multipletests(pvals=se_matrix_2019['p_value'], alpha=0.05, method='bonferroni')[0:2]\n",
    "    se_matrix_2019_sig = se_matrix_2019[se_matrix_2019['sig']==True]\n",
    "    se_matrix_2019_sig.sort_values(by=['SE_rank','p_corrected'], ascending=[True, True])\n",
    "        \n",
    "    \n",
    "    ######Find the SE without any significant drug:\n",
    "    sig_SE = set(list(se_matrix_2019_sig.SE))\n",
    "\n",
    "    ### save the Data Frame\n",
    "    mark = condition[3:]\n",
    "    pickle.dump(se_matrix_2019_sig,  open('../Data/pandemic/top_SE_drug_' + mark+'.pk', 'wb'))\n",
    "    \n",
    "    \n",
    "    ## Save the list of SE that has >1 drug-SE pairs.\n",
    "    pickle.dump(list(set(sig_SE)),  open('../Data/pandemic/sig_drugSE_' + mark+'.pk', 'wb'))\n",
    "#     print(mark, 'kept SE that has sig drug-SE pair:', len(sig_SE))    \n",
    "    \n",
    "    # save the SE removed by drug comfounding\n",
    "    insig_SE = list(set(SE_list) - sig_SE)\n",
    "    pickle.dump(insig_SE,  open('../Data/pandemic/removed_drugSE_' + mark+'.pk', 'wb'))\n",
    "#     print(mark, 'removed SE due to non drug-SE pair:', len(insig_SE))\n",
    "\n",
    "    print('data saved', condition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Drug significantly associated with AE during pandemic\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load data of SE_uncondition_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 512/512 [07:15<00:00,  1.18it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved uncondition_2019_sig_over\n",
      "load data of SE_uncondition_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 187/187 [02:36<00:00,  1.19it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved uncondition_2019_sig_under\n",
      "load data of SE_male_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 135/135 [00:26<00:00,  5.16it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved male_2019_sig_over\n",
      "load data of SE_male_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 25/25 [00:04<00:00,  6.13it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved male_2019_sig_under\n",
      "load data of SE_female_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 281/281 [01:46<00:00,  2.65it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved female_2019_sig_over\n",
      "load data of SE_female_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 78/78 [00:25<00:00,  3.03it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved female_2019_sig_under\n",
      "load data of SE_young_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 54/54 [00:00<00:00, 68.35it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved young_2019_sig_over\n",
      "load data of SE_young_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 14/14 [00:00<00:00, 74.54it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved young_2019_sig_under\n",
      "load data of SE_adult_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 233/233 [00:38<00:00,  6.04it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved adult_2019_sig_over\n",
      "load data of SE_adult_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 16/16 [00:02<00:00,  6.76it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved adult_2019_sig_under\n",
      "load data of SE_elderly_2019_sig_over\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 98/98 [00:11<00:00,  8.21it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved elderly_2019_sig_over\n",
      "load data of SE_elderly_2019_sig_under\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9/9 [00:01<00:00,  8.62it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data saved elderly_2019_sig_under\n"
     ]
    }
   ],
   "source": [
    "# A: drug and SE co-occure in 2020\n",
    "#B: drug appear but SE not appear, B=n_reports -A\n",
    "#C: not drug but SE occure \n",
    "#D: not drug, not SE. C = n_reports -C\n",
    "\n",
    "\n",
    "condition_list = ['SE_uncondition_2019_sig_over',  'SE_uncondition_2019_sig_under', 'SE_male_2019_sig_over', 'SE_male_2019_sig_under',\n",
    "                 'SE_female_2019_sig_over', 'SE_female_2019_sig_under', \n",
    "                 'SE_young_2019_sig_over', 'SE_young_2019_sig_under', 'SE_adult_2019_sig_over', 'SE_adult_2019_sig_under',\n",
    "                 'SE_elderly_2019_sig_over', 'SE_elderly_2019_sig_under']\n",
    "\n",
    "# for condition in condition_list: \n",
    "for condition in condition_list:\n",
    "    pop = condition.split('_')[1]      \n",
    "#     locals()[condition] = pickle.load(open('../Data/pandemic/results/'+condition+'.pk', 'rb'))   \n",
    "    \n",
    "    \n",
    "    # load the data that we want to check\n",
    "    df = pickle.load(open('../Data/pandemic/top_SE_drug_' + condition[3:]+'.pk', 'rb'))\n",
    "    print('load data of',condition)\n",
    "    \n",
    "    df = df.sort_values('p_corrected', ascending=True)\n",
    "    SE_list = list(df.SE)  # find the all the SE\n",
    "    drug_list = list(df.drug) \n",
    "\n",
    "\n",
    "    # for top_SE in top_SE_list:\n",
    "    yr =  2020 #[2019, 2020]\n",
    "\n",
    "    \"\"\"Find the corresponding Drugs of each Top SE\"\"\"\n",
    "    start, end = str(yr) + '-03-11', str(yr) + '-09-31'\n",
    "    ind = [start<i<end for i in all_pd_US_pro['receipt_date']]\n",
    "    input_df= all_pd_US_pro[ind]\n",
    "    n_report = len(input_df)  \n",
    "    \n",
    "    ### split age_grop\n",
    "    input_df['age'] = pd.to_numeric(input_df['age'], errors='coerce')\n",
    "    bins = [1, 20, 65, max(input_df.age)+1]\n",
    "    age_labels = ['Young', 'Adult','Elderly']\n",
    "    input_df['age_group'] = pd.cut(input_df.age, bins, right = False, labels= age_labels)\n",
    "\n",
    "\n",
    "    \"\"\"Limit the reports into specific population!!!!!!!!!!!\"\"\"\n",
    "    if pop =='male':\n",
    "        input_df = input_df[input_df.gender=='1']\n",
    "    elif pop =='female':\n",
    "        input_df = input_df[input_df.gender=='2']\n",
    "    elif pop =='young':\n",
    "        input_df = input_df[input_df.age_group=='Young']\n",
    "    elif pop =='adult':\n",
    "        input_df = input_df[input_df.age_group=='Adult']\n",
    "    elif pop =='elderly':\n",
    "        input_df = input_df[input_df.age_group=='Elderly']\n",
    "\n",
    "    \n",
    "\n",
    "    \"\"\"initialize the Data frame, count the occurance for each SE-drug pairs \"\"\"\n",
    "    se_matrix = pd.DataFrame() #{'SE': SE_list, 'name':list(SE_dic_df.index), 'medra_ID': list(SE_dic_df['medra_ID'])}\n",
    "\n",
    "\n",
    "                    \n",
    "    name = input_df\n",
    "    A, B, C, D =[], [], [],[]\n",
    "    SE_non, drug_non = [], []\n",
    "    for i in tqdm(range(len(SE_list))): # dive into each SE\n",
    "        # Make sure the drug-se pair and the data source\n",
    "        drug = drug_list[i]\n",
    "        se = SE_list[i]           \n",
    "        \n",
    "        # negative and positive reports       \n",
    "        indx_SE = [drug in j for j in name.drugs]  \n",
    "        positive_reports = name[indx_SE]  # The reports that the drug occurred\n",
    "        if len(positive_reports)==0:  # If this drug not occurs in 2020, remove it.\n",
    "            continue\n",
    "        indx_SE_not = [not h for h in indx_SE]\n",
    "        negative_reports = name[indx_SE_not]  \n",
    "        \n",
    "        # Selct the control group\n",
    "        # Extract the independent variables in positive/negative set\n",
    "#         pp =positive_reports.loc[:, 'qualify':'weight']  # features\n",
    "# 'qualify', 'serious', 'receivedate','receiptdate','age','gender','weight','lastingdays'\n",
    "        pp =positive_reports.loc[:, [ 'qualify', 'serious', 'receivedate','receiptdate','age','gender','weight','lastingdays']] \n",
    "\n",
    "        positive_fea = pp.to_numpy()\n",
    "\n",
    "        nn =negative_reports.loc[:, [ 'qualify', 'serious', 'receivedate','receiptdate','age','gender','weight','lastingdays']]\n",
    "        negative_fea = nn.to_numpy()\n",
    "    \n",
    "        scaler = Normalizer()  # MinMaxScaler()# StandardScaler()  # Use normalizer instead of other normalization methods\n",
    "        positive_fea = scaler.fit_transform(positive_fea)\n",
    "        negative_fea = scaler.fit_transform(negative_fea)\n",
    "        \n",
    "        \n",
    "        \"\"\"Trying Faiss\"\"\"\n",
    "        X = positive_fea.astype('float32')\n",
    "        Y = negative_fea.astype('float32')\n",
    "\n",
    "        index = faiss.IndexFlatIP(X.shape[-1])  # Cosine similarity\n",
    "#         Y = Y.copy(order='C')\n",
    "        index.add(np.ascontiguousarray(Y))                  # add negative features;  np.ascontiguousarray(Y)\n",
    "        k = 10                        # we want to see 10 nearest neighbors\n",
    "        _, I = index.search(np.ascontiguousarray(X), k)     # actual search\n",
    "        index_control =  I.reshape([1, -1])\n",
    "        control_group = negative_reports.iloc[index_control[0]]\n",
    "        test_group = positive_reports  \n",
    "    \n",
    "        \n",
    "        # Calculate A, B, C,D for a certain pair\n",
    "        \n",
    "        a,b,c,d =0, 0,0,0\n",
    "        n_test = len(test_group)\n",
    "        n_control = len(control_group)\n",
    "        \n",
    "        \n",
    "        se_ind_test = [se in j for j in test_group.SE]  \n",
    "        se_ind_control = [se in j for j in control_group.SE]  \n",
    "        a = sum(se_ind_test)\n",
    "        c = sum(se_ind_control)\n",
    "        b = n_test -a\n",
    "        d = n_control - c\n",
    "        A.append(a)\n",
    "        B.append(b)\n",
    "        C.append(c)\n",
    "        D.append(d)\n",
    "        \n",
    "        SE_non.append(se)\n",
    "        drug_non.append(drug)\n",
    "\n",
    "    ## len(drug_list) could be smaller than len(drug_non), because some drug paris were droped for unseen in 2020.\n",
    "\n",
    "    se_matrix['SE'] = SE_non\n",
    "    se_matrix['drug'] = drug_non\n",
    "    se_matrix[str(yr)+'_A'] = A\n",
    "    se_matrix[str(yr)+'_B'] = B\n",
    "    se_matrix[str(yr)+'_C'] = C\n",
    "    se_matrix[str(yr)+'_D'] = D\n",
    "\n",
    "    # insert the the drug name and SE name to the dataframe\n",
    "    se_matrix['SE_name'] = se_matrix.apply(lambda row: list(MedDRA_dic_all[MedDRA_dic_all.PT==row.SE]['PT_name'])[0], axis=1)\n",
    "    se_matrix['drug_name'] = se_matrix.apply(lambda row: list(drug_dic_pd[drug_dic_pd['Drugbank_ID']==row.drug].index)[0], axis=1)\n",
    "\n",
    "    \"\"\"Calculate the Fisher's test, and correct the p-value\"\"\"\n",
    "    # Fisher's test\n",
    "#     se_matrix_2019 = se_matrix\n",
    "    se_matrix[str(yr)+'_ROR'] = se_matrix.apply(lambda row: stats.fisher_exact([[row[str(yr)+'_A'], row[str(yr)+'_B']], [row[str(yr)+'_C'], row[str(yr)+'_D']]])[0], axis = 1)\n",
    "    se_matrix[str(yr)+'_ROR_CI_upper'] = se_matrix.apply(lambda row: CI(row[str(yr)+'_ROR'], row[str(yr)+'_A'], row[str(yr)+'_B'],row[str(yr)+'_C'], row[str(yr)+'_D'])[0], axis = 1)\n",
    "    se_matrix[str(yr)+'_ROR_CI_lower'] = se_matrix.apply(lambda row: CI(row[str(yr)+'_ROR'], row[str(yr)+'_A'], row[str(yr)+'_B'],row[str(yr)+'_C'], row[str(yr)+'_D'])[1], axis = 1)\n",
    "    \n",
    "    se_matrix['p_value'] = se_matrix.apply(lambda row: stats.fisher_exact([[row[str(yr)+'_A'], row[str(yr)+'_B']], [row[str(yr)+'_C'], row[str(yr)+'_D']]])[1], axis = 1)\n",
    "\n",
    "    # multipletests\n",
    "    se_matrix['sig'], se_matrix['p_corrected']  = multipletests(pvals=se_matrix['p_value'], alpha=0.05, method='bonferroni')[0:2]\n",
    "    \n",
    "    se_matrix_2019_sig = se_matrix[se_matrix['sig']==True]\n",
    "    \n",
    "    \n",
    "    ### save the Data Frame\n",
    "    mark = condition[3:]\n",
    "    pickle.dump(se_matrix_2019_sig,  open('../Data/pandemic/sig_SE_drug_pair' + mark+ str(yr)+'.pk', 'wb'))\n",
    "    print('data saved', mark)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AEs that satisfy both significance check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load drug-AE pairs during pandemic: uncondition_2019_sig_over\n",
      "load pairs significantly associated with pandemic : uncondition_2019_sig_over\n",
      "kept SE that satisfy both significance check: uncondition_2019_sig_over\n",
      "load drug-AE pairs during pandemic: uncondition_2019_sig_under\n",
      "load pairs significantly associated with pandemic : uncondition_2019_sig_under\n",
      "kept SE that satisfy both significance check: uncondition_2019_sig_under\n",
      "load drug-AE pairs during pandemic: male_2019_sig_over\n",
      "load pairs significantly associated with pandemic : male_2019_sig_over\n",
      "kept SE that satisfy both significance check: male_2019_sig_over\n",
      "load drug-AE pairs during pandemic: male_2019_sig_under\n",
      "load pairs significantly associated with pandemic : male_2019_sig_under\n",
      "kept SE that satisfy both significance check: male_2019_sig_under\n",
      "load drug-AE pairs during pandemic: female_2019_sig_over\n",
      "load pairs significantly associated with pandemic : female_2019_sig_over\n",
      "kept SE that satisfy both significance check: female_2019_sig_over\n",
      "load drug-AE pairs during pandemic: female_2019_sig_under\n",
      "load pairs significantly associated with pandemic : female_2019_sig_under\n",
      "kept SE that satisfy both significance check: female_2019_sig_under\n",
      "load drug-AE pairs during pandemic: young_2019_sig_over\n",
      "load pairs significantly associated with pandemic : young_2019_sig_over\n",
      "kept SE that satisfy both significance check: young_2019_sig_over\n",
      "load drug-AE pairs during pandemic: young_2019_sig_under\n",
      "load pairs significantly associated with pandemic : young_2019_sig_under\n",
      "kept SE that satisfy both significance check: young_2019_sig_under\n",
      "load drug-AE pairs during pandemic: adult_2019_sig_over\n",
      "load pairs significantly associated with pandemic : adult_2019_sig_over\n",
      "kept SE that satisfy both significance check: adult_2019_sig_over\n",
      "load drug-AE pairs during pandemic: adult_2019_sig_under\n",
      "load pairs significantly associated with pandemic : adult_2019_sig_under\n",
      "kept SE that satisfy both significance check: adult_2019_sig_under\n",
      "load drug-AE pairs during pandemic: elderly_2019_sig_over\n",
      "load pairs significantly associated with pandemic : elderly_2019_sig_over\n",
      "kept SE that satisfy both significance check: elderly_2019_sig_over\n",
      "load drug-AE pairs during pandemic: elderly_2019_sig_under\n",
      "load pairs significantly associated with pandemic : elderly_2019_sig_under\n",
      "kept SE that satisfy both significance check: elderly_2019_sig_under\n"
     ]
    }
   ],
   "source": [
    "condition_list = ['SE_uncondition_2019_sig_over',\n",
    "                  'SE_uncondition_2019_sig_under', 'SE_male_2019_sig_over', 'SE_male_2019_sig_under',\n",
    "                 'SE_female_2019_sig_over', 'SE_female_2019_sig_under', \n",
    "                 'SE_young_2019_sig_over', 'SE_young_2019_sig_under', 'SE_adult_2019_sig_over', 'SE_adult_2019_sig_under',\n",
    "                 'SE_elderly_2019_sig_over', 'SE_elderly_2019_sig_under']\n",
    "\n",
    "for condition in condition_list: \n",
    "    mark = condition[3:]\n",
    "    yr = 2020    \n",
    "    yr_pair = pickle.load(open('../Data/pandemic/sig_SE_drug_pair' + mark+ str(yr)+'.pk', 'rb'))\n",
    "    print('load drug-AE pairs during pandemic: {}'.format(mark))\n",
    "    \n",
    "    ### Load the pandemic-sig pairs\n",
    "    pandemic_pair = pickle.load(open('../Data/pandemic/top_SE_drug_' + mark+'.pk', 'rb'))    \n",
    "    print('load pairs significantly associated with pandemic : {}'.format(mark))\n",
    "\n",
    "    \n",
    "    ## Merge them  to get both_sig_pair\n",
    "    both_sig_pair = pd.merge(pandemic_pair, yr_pair, how='inner', on=['SE','drug'])\n",
    "    \n",
    "    ## Save the list of SE that has >1 drug-SE pairs.\n",
    "    sig_SE = list(both_sig_pair.SE)\n",
    "    pickle.dump(list(set(sig_SE)),  open('../Data/pandemic/sig_drugSE_' + mark+'.pk', 'wb'))\n",
    "    print( 'kept SE that satisfy both significance check:', mark)   "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "274.358px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}