manuscript_phenotype_script.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:09:54.848379Z",
     "start_time": "2022-01-13T21:09:53.504204Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import csv\n",
    "import paramiko\n",
    "import gzip\n",
    "\n",
    "import pymysql\n",
    "from sshtunnel import SSHTunnelForwarder\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "from datetime import timedelta, date, datetime\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "from collections import defaultdict\n",
    "import re\n",
    "from shutil import copyfile\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:00.870777Z",
     "start_time": "2022-01-13T21:10:00.858969Z"
    }
   },
   "outputs": [],
   "source": [
    "#Inputs\n",
    "\n",
    "database = \"clinical_merge_v5_2020q4\"\n",
    "new_rels = True\n",
    "\n",
    "#Build from ccsr\n",
    "ped_data = '/data1/home/kel2158/Projects/pedigree_v3/data/final_results/'\n",
    "output_folder= '/home/kel2158/Projects/COVIDHerit/data/12_29_21/'\n",
    "hhid_folder = 'home/kel2158/Projects/hhid_update/data/'\n",
    "\n",
    "#Input files\n",
    "finename_fam_input = ped_data +\"family_v3.txt.gz\"\n",
    "filename_hhid =  hhid_folder+'hhid_trimmed_notcln.txt.gz'\n",
    "\n",
    "#Output filenames\n",
    "filename_demo = output_folder+ 'demo.txt'\n",
    "filename_covariates = output_folder+'covariates_building_race.txt'\n",
    "filename_covariates_raceonly = output_folder+'covariates_race.txt'\n",
    "\n",
    "filename_pheno_binary = output_folder+'covid_binary.txt'\n",
    "filename_pheno_quant = output_folder+'covid_quant.txt'\n",
    "\n",
    "\n",
    "fam_count_file = output_folder + \"/fam_size.txt\"\n",
    "\n",
    "#Variables\n",
    "#Patient must have died within 90 days of latest positve test\n",
    "death_cutoff = datetime.timedelta(days=90)\n",
    "intub_cutoff = datetime.timedelta(days=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:02.491867Z",
     "start_time": "2022-01-13T21:10:02.482026Z"
    }
   },
   "outputs": [],
   "source": [
    "def return_sign_in(login_loc = '/home/kel2158/login.txt'):\n",
    "    login_dict = dict()\n",
    "    with open(login_loc, \"r\") as file:\n",
    "        for line in file.readlines():\n",
    "            login_dict[ line.split(\",\")[0]] = line.split(\",\")[1].strip()\n",
    "    return(login_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:03.017935Z",
     "start_time": "2022-01-13T21:10:02.960350Z"
    }
   },
   "outputs": [],
   "source": [
    "login_dict = return_sign_in()\n",
    "\n",
    "ssh = paramiko.SSHClient()\n",
    "ssh.set_missing_host_key_policy(paramiko.AutoAddPolicy())\n",
    "ssh.connect('10.144.220.21', username=login_dict['user'], password=login_dict['mimir_pw'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:04.657807Z",
     "start_time": "2022-01-13T21:10:04.443279Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Make sure to call server.stop() when you want to disconnect\n",
    "# after calling db.close()\n",
    "server = SSHTunnelForwarder('10.144.220.21',\n",
    "                            ssh_username='kel2158',\n",
    "                            ssh_password=login_dict['mimir_pw'],\n",
    "                            remote_bind_address=('127.0.0.1', 3306))\n",
    "server.start()\n",
    "\n",
    "db = pymysql.connect(host='127.0.0.1',\n",
    "                     port=server.local_bind_port,\n",
    "                     user=login_dict['user'],\n",
    "                     password=login_dict[\"db_pw\"],\n",
    "                     db=database)\n",
    "\n",
    "cur = db.cursor()\n",
    "\n",
    "cur.execute( \"SET sql_mode=(SELECT REPLACE(@@sql_mode,'ONLY_FULL_GROUP_BY',''));\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull lab and family data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:12.381741Z",
     "start_time": "2022-01-13T21:10:08.821426Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in family file \n",
    "file_name = finename_fam_input\n",
    "ftp = ssh.open_sftp()\n",
    "file = ftp.file(file_name, \"r+\")\n",
    "gfile = gzip.open(file)\n",
    "\n",
    "fam_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "ftp.close()\n",
    "ssh.close()\n",
    "mrn2famid = dict(zip (fam_df[\"ptid\"],fam_df[\"famid\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-29T21:47:24.084191Z",
     "start_time": "2021-12-29T21:47:01.244853Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "496876"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Pull patients who have had positive covid tests\n",
    "\n",
    "pull_data_query = \"\"\"\n",
    "DROP TABLE IF EXISTS user_kel2158.covid_lab_table;\n",
    "\"\"\"\n",
    "cur.execute(pull_data_query)\n",
    "\n",
    "\n",
    "pull_data_query = \"\"\" CREATE TABLE user_kel2158.covid_lab_table AS\n",
    "\tSELECT  DISTINCT l.pat_mrn_id, m.order_time, l.order_proc_id,  l.result_date, l.ord_value, l.person_id_2020, l.person_id_2019,\n",
    "\tCASE\n",
    "\tWHEN l.ord_value IN (\"Detected\",\n",
    "        \"Detected A\",\n",
    "        \"Detected..\",\n",
    "        \"POSITIVE\",\n",
    "        \"Postive\",\n",
    "        \"Postivie A\",\n",
    "        \"Yes\"\n",
    "        ) THEN 'pos'\n",
    "\tWHEN l.ord_value IN (\"Not Detected\",\n",
    "        \"Negative\",\n",
    "        \"No\",\n",
    "        \"Not Detected (External)\",\n",
    "        \"SARS-CoV-2 neg\",\n",
    "        \"neg\",\n",
    "        \"Non detected\",\n",
    "        \"Undetected\",\n",
    "        \"Not Dectected\",\n",
    "        \"NotDetectd\",\n",
    "        \"Negatvie\",\n",
    "        \"NOT_DETECT\",\n",
    "        \"notdetected\",\n",
    "        \"Not Deteceted\",\n",
    "        \"Note Detected\",\n",
    "        \"not dtected\",\n",
    "        \"not  detected\") THEN 'neg'\n",
    "\tEND AS test_result\n",
    "    FROM covid_labs l,\n",
    "    covid_measurements m\n",
    "    WHERE \n",
    "\t(m.procedure_name like '%SARS%'\n",
    "\tOR m.component_name like '%SARS%'\n",
    "    OR m.component_loinc_code = 94500)  AND m.order_proc_id = l.order_proc_id\n",
    ";\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "\n",
    "pull_data_query = \"CREATE INDEX pat_mrn_id ON user_kel2158.covid_lab_table (pat_mrn_id);\"\n",
    "cur.execute(pull_data_query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:13.097421Z",
     "start_time": "2022-01-13T21:10:12.816909Z"
    }
   },
   "outputs": [],
   "source": [
    "#Get patients with an icd code for history of covid \n",
    "pull_data_query = \"\"\"\n",
    "SELECT DISTINCT p.pat_mrn_id\n",
    "FROM clinical_merge_v5_2020q4.covid_patients p\n",
    "WHERE icd10_code = 'Z86.16';\n",
    "\n",
    "\"\"\"\n",
    "cur.execute(pull_data_query)\n",
    "\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "hist_cov_ptid = set([e[0] for e in rawdata])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:15.649538Z",
     "start_time": "2022-01-13T21:10:13.968282Z"
    }
   },
   "outputs": [],
   "source": [
    "#Pull patients who have had positive covid tests\n",
    "\n",
    "pull_data_query = \"\"\"SELECT DISTINCT pat_mrn_id, order_time\n",
    "FROM user_kel2158.covid_lab_table\n",
    "WHERE test_result = 'pos';\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "p_df = pd.DataFrame(list( rawdata ),columns=[\"ptid\",\"order_time\"])\n",
    "\n",
    "cases = set([e[0] for e in rawdata])\n",
    "\n",
    "p_df[\"value\"] = 1\n",
    "\n",
    "p_df[\"famid\"] = p_df[\"ptid\"].map(mrn2famid)\n",
    "\n",
    "p_df = p_df[ p_df[\"famid\"].notnull() ] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:19.624636Z",
     "start_time": "2022-01-13T21:10:16.251315Z"
    }
   },
   "outputs": [],
   "source": [
    "#Pull patients who have had negative covid tests\n",
    "\n",
    "pull_data_query = \"\"\"SELECT DISTINCT pat_mrn_id\n",
    "FROM user_kel2158.covid_lab_table\n",
    "WHERE test_result = 'neg'\n",
    ";\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "intermediate_controls = set([e[0] for e in rawdata])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:10:46.004286Z",
     "start_time": "2022-01-13T21:10:36.918730Z"
    }
   },
   "outputs": [],
   "source": [
    "#Get familial degrees\n",
    "pull_data_query = \"\"\"\n",
    "SELECT r.mrn, r.relation_mrn, d.relationship, d.degree\n",
    "FROM clinical_relationships_v3.all_relationships_to_generate_pedigree_file r,\n",
    "user_kel2158.covid_lab_table c1,\n",
    "user_kel2158.covid_lab_table c2,\n",
    "clinical_relationships.relationships_degree d\n",
    "WHERE r.mrn = c1.pat_mrn_id AND\n",
    "r.relation_mrn = c2.pat_mrn_id AND\n",
    "r.relationship = d.relationship\n",
    ";\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "relationship_df = pd.DataFrame(rawdata,columns=[\"mrn\",\"rel_mrn\",\"relationship\",\"degree\"])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:19.669652Z",
     "start_time": "2022-01-13T21:11:19.658145Z"
    }
   },
   "outputs": [],
   "source": [
    "def pull_fam_count( pheno, df, print_res = True, \\\n",
    "                   return_fams_df=False, type_pheno=\"binary_withproband\",\n",
    "                   return_patient_limited_df = False ):\n",
    "    \n",
    "    \"\"\"Returns the number of families in a phenotype for SOLAR based on pheno \n",
    "     Needs phenotype column, famid column and ptid column\"\"\"    \n",
    "    \n",
    "    #subset on phenotype of interest and those in our pedigree\n",
    "    df_phenosub = df[ (df[\"pheno\"] == pheno ) \n",
    "                       &  (~df[\"famid\"].isnull())] \n",
    "    \n",
    "    #Get the total number of  members in each fam and then number of cases\n",
    "    df_phenogrouped = df_phenosub.groupby( 'famid' )\\\n",
    "        .agg({ \"ptid\":\"count\",\"value\":sum })\\\n",
    "        .reset_index()\\\n",
    "        .rename(columns={ \"ptid\":\"num_fam_members\",\"value\":\"num_cases\" })\n",
    "    \n",
    "    #Subset to fams that will be included in SOLAR\n",
    "    #ie has more than one family member an at least one case in each family\n",
    "    if type_pheno == \"binary_withproband\":\n",
    "        solar_fams_df = df_phenogrouped[ (df_phenogrouped[\"num_fam_members\"] > 1) \\\n",
    "                                        & (df_phenogrouped[\"num_cases\"] > 0) ]\n",
    "    \n",
    "    else:\n",
    "         solar_fams_df = df_phenogrouped[ (df_phenogrouped[\"num_fam_members\"] > 1)]\n",
    "    \n",
    "    solar_limited_df = df_phenosub[ df_phenosub[\"famid\"].isin(solar_fams_df[\"famid\"].unique())]\n",
    "    solar_limited_df.drop_duplicates(inplace=True)\n",
    "    solar_patients = set( solar_limited_df[\"ptid\"] )\n",
    "    \n",
    "    \n",
    "    \n",
    "    if print_res:\n",
    "            print(f'There are {len(solar_fams_df)} SOLAR families')\n",
    "    \n",
    "    if return_fams_df:\n",
    "        return solar_fams_df, solar_patients\n",
    "    \n",
    "    elif return_patient_limited_df == True:\n",
    "        return solar_limited_df\n",
    "    \n",
    "    else:\n",
    "        return(len( solar_fams_df) )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:24.965965Z",
     "start_time": "2022-01-13T21:11:24.953873Z"
    }
   },
   "outputs": [],
   "source": [
    "def write_famcount(count, pheno):\n",
    "    \n",
    "    orig_file = fam_count_file\n",
    "    temp_file = fam_count_file[:-4]+\"_temp.txt\"\n",
    "    \n",
    "    #Copy the original file into a temp file\n",
    "    copyfile(orig_file, temp_file)  \n",
    "        \n",
    "    #Delete original file\n",
    "    os.remove(orig_file)\n",
    "    \n",
    "    bad_words = [pheno]\n",
    "    \n",
    "    #exclude pheno if it's already there\n",
    "    with open(temp_file) as oldfile, open(orig_file, 'w') as newfile:\n",
    "        for line in oldfile:\n",
    "            if not any(bad_word in line for bad_word in bad_words):\n",
    "                newfile.write(line)\n",
    "    \n",
    "    #append pheno to end of updated file\n",
    "    with open(orig_file, \"a\") as file_object:\n",
    "        #append pheno to end of file\n",
    "        file_object.write(pheno+\" \"+str(count))\n",
    "        file_object.write(\"\\n\")\n",
    "    \n",
    "    file_object = open(orig_file, 'a')\n",
    "    \n",
    "    #Delete temporary file\n",
    "    os.remove(temp_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheno: Infection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:28.786584Z",
     "start_time": "2022-01-13T21:11:27.634226Z"
    }
   },
   "outputs": [],
   "source": [
    "#Controls Overall: had a test but never positive\n",
    "controls = intermediate_controls-cases\n",
    "\n",
    "#Remove those patients with history of covid icd code\n",
    "#Note, don't use covid only icd codes because gold standard is pcr and need timing of infection\n",
    "controls = controls - hist_cov_ptid\n",
    "controls = set(controls)\n",
    "\n",
    "pheno_df = pd.concat( [  pd.DataFrame(zip(cases,len(cases)*[\"covid\"],len(cases)*[1])),\n",
    "                       pd.DataFrame(zip(controls,len(controls)*['covid'],len(controls)*[0])) ])\n",
    "pheno_df.rename(columns={0:\"ptid\",1:\"pheno\",2:\"value\"},inplace=True)\n",
    "\n",
    "#Drop duplicates not necessary, select distinct negatives and select min and max for positives\n",
    "\n",
    "\n",
    "pheno_df[\"famid\"] = pheno_df[\"ptid\"].map(mrn2famid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:31.560600Z",
     "start_time": "2022-01-13T21:11:31.552367Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "181397"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(controls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T20:55:50.685469Z",
     "start_time": "2022-01-13T20:55:50.677640Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17875"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(cases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:51.697857Z",
     "start_time": "2022-01-13T21:11:37.145349Z"
    }
   },
   "outputs": [],
   "source": [
    "pull_data_query = \"\"\" SELECT DISTINCT pat_mrn_id, death_date\n",
    "FROM covid_persons\n",
    "WHERE death_date IS NULL; \"\"\"\n",
    "\n",
    "\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "deaths_df = pd.DataFrame(list(rawdata),columns=[\"mrn\",\"death_date\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T03:48:44.292156Z",
     "start_time": "2021-12-30T03:48:44.142508Z"
    }
   },
   "outputs": [],
   "source": [
    "#Output file file\n",
    "file_name = output_folder+'covid.txt'\n",
    "\n",
    "pheno_df[ pheno_df[\"famid\"].notnull() ][cols_to_keep].drop_duplicates().to_csv\\\n",
    "(file_name,index=False, sep='\\t')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T16:51:40.466451Z",
     "start_time": "2021-12-30T16:51:40.410019Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 1324 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"covid\", pheno_df),\\\n",
    "              \"covid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T03:48:58.378935Z",
     "start_time": "2021-12-30T03:48:58.205225Z"
    }
   },
   "outputs": [],
   "source": [
    "pheno_df[\"pheno\"] = \"covid_noproband\"\n",
    "\n",
    "file_name = output_folder+'covid_noproband.txt'\n",
    "pheno_df[ pheno_df[\"famid\"].notnull() ][cols_to_keep].drop_duplicates().to_csv\\\n",
    "(file_name,index=False, sep='\\t')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-18T20:07:01.128667Z",
     "start_time": "2022-01-18T20:07:01.060567Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 5676 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "pheno_df[\"pheno\"] = \"covid_noproband\"\n",
    "write_famcount( pull_fam_count(\"covid_noproband\", pheno_df, type_pheno=\"noproband\"),\\\n",
    "              \"covid_noproband\")\n",
    "\n",
    "pheno_df[\"pheno\"] = \"covid\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High Degree Families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:53.745794Z",
     "start_time": "2022-01-13T21:11:53.333289Z"
    }
   },
   "outputs": [],
   "source": [
    "high_deg_fams = pheno_df[ pheno_df[\"ptid\"].astype(str).isin\\\n",
    "         (set( relationship_df[ relationship_df[\"degree\"].isin\\\n",
    "                               ([\"Second\",\"Third\",\"Fourth\"]) ][\"mrn\"]) )][\"famid\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:11:55.456045Z",
     "start_time": "2022-01-13T21:11:55.441734Z"
    }
   },
   "outputs": [],
   "source": [
    "#Limit to higher degree\n",
    "covid_notfirstdegree_df = pheno_df[ pheno_df[\"famid\"].isin(high_deg_fams) ]\n",
    "\n",
    "covid_notfirstdegree_df[\"pheno\"] = \"covid_highdeg_fam\"\n",
    "\n",
    "#Create no proband requirement\n",
    "covid_notfirstdegree_nopro_df = covid_notfirstdegree_df.copy()\n",
    "covid_notfirstdegree_nopro_df[\"pheno\"] = \"covid_highdeg_nopro\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-18T20:08:29.918628Z",
     "start_time": "2022-01-18T20:08:29.898486Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 296 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"covid_highdeg_fam\", covid_notfirstdegree_df),\\\n",
    "              \"covid_highdeg_fam\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-18T20:07:55.315611Z",
     "start_time": "2022-01-18T20:07:55.291156Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 962 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"covid_highdeg_nopro\", covid_notfirstdegree_nopro_df, type_pheno=\"noproband\"),\\\n",
    "               \"covid_highdeg_nopro\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T03:49:11.882368Z",
     "start_time": "2021-12-30T03:49:11.842482Z"
    }
   },
   "outputs": [],
   "source": [
    "#Output file file\n",
    "file_name = output_folder+'covid_high_degfam.txt'\n",
    "covid_notfirstdegree_df[cols_to_keep].drop_duplicates().to_csv\\\n",
    "(file_name,index=False, sep='\\t')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T03:49:18.905458Z",
     "start_time": "2021-12-30T03:49:18.881704Z"
    }
   },
   "outputs": [],
   "source": [
    "file_name = output_folder+'covid_high_degfam_noproband.txt'\n",
    "covid_notfirstdegree_nopro_df[cols_to_keep].drop_duplicates().to_csv(file_name,index=False, sep='\\t')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull Time Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:49:17.304811Z",
     "start_time": "2022-01-13T21:49:15.994493Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03731468531468531"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Get earliest PCR! \n",
    "pull_data_query = \"\"\" SELECT DISTINCT pat_mrn_id, MIN(order_time) AS order_time, MAX(order_time)\n",
    "                      FROM user_kel2158.covid_lab_table\n",
    "                      WHERE test_result = \"pos\" \n",
    "                    GROUP BY pat_mrn_id;  \"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "first_pcr_nohospreq_df = pd.DataFrame(list( rawdata ), columns=[\"mrn\",\n",
    "                           \"order_time\",\"latest_pcr_order_time\"\n",
    "                           ] )\n",
    "\n",
    "#Get time difference between first and last positive test\n",
    "first_pcr_nohospreq_df[\"pcr_timediff\"] = (first_pcr_nohospreq_df[\"latest_pcr_order_time\"] \\\n",
    "                         - first_pcr_nohospreq_df[\"order_time\"] ).apply( lambda x: x.days ) \n",
    "\n",
    "\n",
    "len( first_pcr_nohospreq_df[ first_pcr_nohospreq_df[\"pcr_timediff\"] > 30 ] )\n",
    "\n",
    "#Percent patients with a second covid test after 30 days...\n",
    "len( first_pcr_nohospreq_df[ first_pcr_nohospreq_df[\"pcr_timediff\"] > 30 ] ) / len( first_pcr_nohospreq_df )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:49:20.600889Z",
     "start_time": "2022-01-13T21:49:20.590868Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "#Every case should have an earlier positive amount\n",
    "len( first_pcr_nohospreq_df ) == len(cases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:49:24.664709Z",
     "start_time": "2022-01-13T21:49:23.753603Z"
    }
   },
   "outputs": [],
   "source": [
    "#Order by eariliest test\n",
    "first_pcr_nohospreq_df.sort_values(by=\"order_time\", inplace=True)\n",
    "\n",
    "#A control is same at any time cut (neg test and no positive ever)\n",
    "#Set pheno value 1 for cases (infected), 0 for controls (never infected)\n",
    "first_pcr_nohospreq_df[\"value\"] = 1\n",
    "\n",
    "#Turn datetime into date\n",
    "first_pcr_nohospreq_df[\"order_date\"] = first_pcr_nohospreq_df[\"order_time\"].apply(lambda x: x.date())\n",
    "\n",
    "#Map famids to mrns\n",
    "first_pcr_nohospreq_df[\"famid\"] = first_pcr_nohospreq_df[\"mrn\"].map( mrn2famid )\n",
    "\n",
    "#Rename columns to align with pull family function\n",
    "first_pcr_nohospreq_df.rename(columns={\"mrn\":\"ptid\"},inplace=True)\n",
    "first_pcr_nohospreq_df.drop_duplicates(inplace=True)\n",
    "#Remove those that aren't in fam file\n",
    "first_pcr_nohospreq_df = first_pcr_nohospreq_df[ ~first_pcr_nohospreq_df[\"famid\"].isnull() ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull Admissions Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:51:34.795377Z",
     "start_time": "2022-01-13T21:50:54.263292Z"
    }
   },
   "outputs": [],
   "source": [
    "#Get first admission that is within 96 hours of test\n",
    "#And \n",
    "pull_data_query = \"\"\"\n",
    "\n",
    "SELECT *\n",
    "FROM (\n",
    "    SELECT c.pat_mrn_id, a.encounter_start_datetime, a.encounter_censor_datetime,\n",
    "    EXTRACT(HOUR FROM TIMEDIFF(l.order_time, a.encounter_start_datetime) ) AS hours_from_admission, \n",
    "    l.order_time, a.pat_enc_csn_id, c.event_type_desc\n",
    "\n",
    "    FROM covid_ADT_defined_visits a,\n",
    "\tcovid_ADT c,\n",
    "        \n",
    "        #Select first positive covid tests that's not an antigen test!\n",
    "        ( SELECT DISTINCT pat_mrn_id, MIN(order_time) AS order_time\n",
    "                      FROM user_kel2158.covid_lab_table\n",
    "                      WHERE test_result = \"pos\" \n",
    "                    GROUP BY pat_mrn_id ) l\n",
    "\n",
    "  WHERE c.pat_mrn_id = l.pat_mrn_id \n",
    "  AND c.pat_enc_csn_id = a.pat_enc_csn_id) d\n",
    "\n",
    "WHERE\n",
    "    #Covid test is within  96 hours of \"admission\", doesn't matter how long patient in hosp for, can be outpatient\n",
    "    0 <= d.hours_from_admission \n",
    "    AND 96 >= d.hours_from_admission\n",
    ";\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "first_admsn_df = pd.DataFrame(list( rawdata ), columns=[\"mrn\",\n",
    "                           \"admsn_time\",\"disch_time\",\"hrs_from_admission\"\n",
    "                            ,\"first_pos_result_datetime\", \"pat_enc\",\"visit_type\"\n",
    "                           ] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:51:55.899456Z",
     "start_time": "2022-01-13T21:51:38.355504Z"
    }
   },
   "outputs": [],
   "source": [
    "#Pull all admissions for patients with a positive pcr test\n",
    "pull_data_query = \"\"\"\n",
    "               SELECT DISTINCT c.pat_mrn_id, a.encounter_start_datetime, a.encounter_censor_datetime, \n",
    "                l.order_time, a.pat_enc_csn_id, a.encounter_censor_type\n",
    "\n",
    "                FROM covid_ADT_defined_visits a,\n",
    "                covid_ADT c,\n",
    "\n",
    "        #Select first positive covid tests that's not an antigen test!\n",
    "        ( SELECT DISTINCT pat_mrn_id, MIN(order_time) AS order_time\n",
    "                      FROM user_kel2158.covid_lab_table\n",
    "                      WHERE test_result = \"pos\" \n",
    "                    GROUP BY pat_mrn_id ) l\n",
    "\n",
    "              WHERE c.pat_mrn_id = l.pat_mrn_id\n",
    "              AND c.pat_enc_csn_id = a.pat_enc_csn_id\n",
    "             \n",
    ";\"\"\"\n",
    "\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "admsn_df = pd.DataFrame(list( rawdata ),columns = [\"mrn\",\"hosp_admsn\",\n",
    "                                                 \"hosp_disch\",\"first_pos_result_datetime\",\n",
    "                                                 \"pat_enc_csn_id\",\"censor_stat\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:51:58.221162Z",
     "start_time": "2022-01-13T21:51:58.198760Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Discharge', 'Hospital Outpatient', 'Census',\n",
       "       'Leave of Absence Census'], dtype=object)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "admsn_df[\"censor_stat\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:52:00.511959Z",
     "start_time": "2022-01-13T21:52:00.474582Z"
    }
   },
   "outputs": [],
   "source": [
    "#Closed admissions\n",
    "discharged_admissions = admsn_df[ admsn_df[\"censor_stat\"].isin\\\n",
    "                             (['Discharge', 'Hospital Outpatient']) ]\n",
    "\n",
    "#Note admissions that are still open, will exclude for patients w/ days hospitalized potentially\n",
    "open_admission = admsn_df[ admsn_df[\"censor_stat\"].isin( \\\n",
    "                 [\"Census\", 'Leave of Absence Census' ])\n",
    "                          & \n",
    "                         ~admsn_df[\"pat_enc_csn_id\"].isin(discharged_admissions[\"pat_enc_csn_id\"]) \n",
    "                         ]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:52:02.703910Z",
     "start_time": "2022-01-13T21:52:02.696592Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "161"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(open_admission)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pull DNI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:26:37.819046Z",
     "start_time": "2021-12-30T17:12:29.773549Z"
    }
   },
   "outputs": [],
   "source": [
    "# Individuals with DNI orders\n",
    "pull_data_query =\"\"\"\n",
    "SELECT pat_mrn_id, MIN(order_date) AS dnr_date\n",
    "FROM covid_orders\n",
    "WHERE description = 'DNR/DNI-DO NOT RESUSCITATE/DO NOT INTUBATE'\n",
    "GROUP BY pat_mrn_id;\n",
    " \"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "mrn2dni_dnr = {e[0]:e[1] for e in rawdata}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:26:54.147559Z",
     "start_time": "2021-12-30T17:26:52.612604Z"
    }
   },
   "outputs": [],
   "source": [
    "#Pull intubation data\n",
    "pull_data_query = \"\"\"\n",
    "SELECT DISTINCT pat_mrn_id,order_date\n",
    "FROM covid_intubation_orders\n",
    "WHERE order_date IS NOT NULL AND order_status != 'Canceled';\"\"\"\n",
    "     \n",
    "cur.execute(pull_data_query)\n",
    "\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "intub_df = pd.DataFrame(list(rawdata),columns=[\"mrn\",\"intubation_date\"] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo and Covariates\n",
    "### Pull Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:30:30.425659Z",
     "start_time": "2021-12-30T17:27:09.410532Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "647231\n",
      "438774\n",
      "647231\n"
     ]
    }
   ],
   "source": [
    "#Pull addresses from covid\n",
    "\n",
    "pull_data_query = \"\"\"\n",
    " SELECT pat_mrn_id, add_line_1,zip\n",
    "FROM covid_persons\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "print ( len(rawdata) )\n",
    "add_df_covid = pd.DataFrame(list(rawdata),\n",
    "                columns=[\"mrn\",\"address\",\"zip\"])\n",
    "\n",
    "\n",
    "#Update to only include patients we have lab covid data for, since hhid file already created\n",
    "#Note, no cumc_demogs data with address in 2020q4 db\n",
    "pull_data_query = \"\"\"\n",
    " SELECT mrn, add_line_1,zip\n",
    "FROM clinical_merge_v5_240919.cumc_patient_demogs\n",
    "WHERE mrn IN (\n",
    "    SELECT DISTINCT l.pat_mrn_id\n",
    "    FROM covid_labs l)\n",
    "\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "print ( len(rawdata) )\n",
    "\n",
    "add_df_2019 = pd.DataFrame(list(rawdata),\n",
    "                columns=[\"mrn\",\"address\",\"zip\"])\n",
    "\n",
    "\n",
    "# First try using covid data (more up to date!)\n",
    "add_df = pd.concat([ add_df_2019[ ~add_df_2019[\"mrn\"].isin(add_df_covid[\"mrn\"]) ], add_df_covid ])\n",
    "\n",
    "add_df[\"address\"].update( add_df[\"address\"].str.upper() )\n",
    "\n",
    "add_df.drop_duplicates(inplace=True)\n",
    "\n",
    "print(len(add_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:31:44.067546Z",
     "start_time": "2021-12-30T17:30:48.773633Z"
    }
   },
   "outputs": [],
   "source": [
    "#1 codes for \"Female\"\n",
    "#Pull demo data from persons table if missing\n",
    "\n",
    "\n",
    "pull_data_query = \"\"\"SELECT p.pat_mrn_id, p.sex_c, p.race_1, p.ethnicity, YEAR(p.birth_date)\n",
    "FROM covid_persons p,\n",
    "    covid_labs l\n",
    "WHERE p.pat_mrn_id = l.pat_mrn_id;\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "\n",
    "pull_data_query = \"\"\"SELECT p2.pat_mrn_id, gender_source_value, race_source_value, ethnicity_source_value, year_of_birth\n",
    "FROM person p,\n",
    "    covid_patient2person p2,\n",
    "    covid_labs l\n",
    "WHERE \n",
    "    l.pat_mrn_id = p2.pat_mrn_id\n",
    "    AND p2.person_id_2020 = p.person_id\"\"\"\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata1 = cur.fetchall()\n",
    "\n",
    "\n",
    "\n",
    "#Dedup data\n",
    "rawdata = set(rawdata)\n",
    "rawdata1 = set(rawdata1)\n",
    "\n",
    "#Add in covid demo data first\n",
    "demo_covid_df = pd.DataFrame(list(rawdata), columns=[\"mrn\",\"sex_raw\",\"race1\",\"ethnicity\",\"birth_year\"])\n",
    "\n",
    "#Supplement with cumc demog data \n",
    "\n",
    "supplement_demo_df = pd.DataFrame(list(rawdata1), columns=[\"mrn\",\"sex_raw\",\"race1\",\"ethnicity\",\"birth_year\"])\n",
    "\n",
    "demo_df = pd.concat([supplement_demo_df[ ~supplement_demo_df['mrn'].\\\n",
    "                                        isin(demo_covid_df[\"mrn\"].unique())] ,\n",
    "                   demo_covid_df])\n",
    "\n",
    "demo_df.drop_duplicates(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 327,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:35:04.696252Z",
     "start_time": "2021-12-30T17:35:04.301222Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "set()"
      ]
     },
     "execution_count": 327,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Are we adding anything with supplementary data?\n",
    "#False if we are adding anything\n",
    "len( set( demo_covid_df['mrn'] ).union( set( supplement_demo_df['mrn'] ) ) ) == len( set( demo_covid_df['mrn'] ))\n",
    "\n",
    "#Are we adding anything with supplementary data?\n",
    "#Empty set if no\n",
    "set( supplement_demo_df['mrn'] ) -  ( set( demo_covid_df['mrn'] ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean Address Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 328,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:36:45.626452Z",
     "start_time": "2021-12-30T17:35:22.967516Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "647231it [01:21, 7906.51it/s]\n"
     ]
    }
   ],
   "source": [
    "#Shares house with someone who tests positive\n",
    "address2mrnv2 = defaultdict(list)\n",
    "mrn2address = dict()\n",
    "\n",
    "add_inall_df =  add_df\n",
    "add_inall_df.drop_duplicates(inplace=True)\n",
    "add_inall_df.fillna(\"(null)\",inplace=True)\n",
    "\n",
    "for  row in tqdm(add_inall_df.iterrows()):\n",
    "    #First row header\n",
    "    if row[0] != 0:\n",
    "        mrn = row[1][\"mrn\"]\n",
    "        #Replace multiple spaces with one space\n",
    "        address = re.sub( r'\\s{2}' , \" \", row[1][\"address\"] ) \n",
    "        house_num = address.split(\" \")[0]\n",
    "        #We want only the first part of zip since not all have full\n",
    "        zipc = row[1][\"zip\"].split(\"-\")[0]\n",
    "        \n",
    "        try:\n",
    "            street =  address.replace(\" W \",\" \")\\\n",
    "                                    .replace(\" W. \",\" \")\\\n",
    "                                    .replace(\" WEST \",\" \")\\\n",
    "                                    .replace(\" E \",\" \")\\\n",
    "                                    .replace(\" E. \",\" \")\\\n",
    "                                    .replace(\" EAST \",\" \").split(\" \")[1]\n",
    "            hid = house_num+\"_\"+street+\"_\"+zipc\n",
    "            \n",
    "            #If missing a core component of hid or po box, don't add\n",
    "            if ( (hid.find(\"null\") == -1) & \n",
    "                (hid.find(\"__\") == -1) &\n",
    "                (hid.find(\"BOX_\") == -1)& \n",
    "                (hid.find(\"PO_BOX\") == -1) &\n",
    "                (hid.find(\"POBOX\") == -1) \n",
    "               ) :\n",
    "                address2mrnv2[hid].append(mrn)\n",
    "                mrn2address[mrn] = hid\n",
    "\n",
    "        #No space to parse out street number\n",
    "        #Updated to exclude, would rather have\n",
    "        except:\n",
    "            #hid = house_num+\"_\"+zipc\n",
    "            continue\n",
    "\n",
    "#         address2mrnv2[hid].append(mrn)\n",
    "#         mrn2address[mrn] = hid\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:37:10.768344Z",
     "start_time": "2021-12-30T17:37:09.762256Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████| 283950/283950 [00:00<00:00, 286733.31it/s]\n"
     ]
    }
   ],
   "source": [
    "#List of buildings with more than one individual in it tested where \n",
    "#at least one person has covid\n",
    "hid2covidcnt = dict()\n",
    "mrn2covidcnt =  dict()\n",
    "\n",
    "#For each building with mrns\n",
    "for hid, mrns in tqdm( address2mrnv2.items() ):\n",
    "    covid_count = 0\n",
    "   \n",
    "    #For each patient in building, check if had covid\n",
    "    for mrn in mrns:\n",
    "        if mrn in cases:\n",
    "            covid_count += 1\n",
    "\n",
    "    hid2covidcnt[hid] = covid_count\n",
    "    \n",
    "    #How many people in building other than self have covid?\n",
    "    for mrn in mrns:\n",
    "        if covid_count == 0:\n",
    "            mrn2covidcnt[mrn] = covid_count\n",
    "        \n",
    "        #if patient has covid, remove one from count and assign \n",
    "        elif mrn in cases:\n",
    "            mrn2covidcnt[mrn] = covid_count-1\n",
    "        \n",
    "        #if patient doesn't have covid assign total count\n",
    "        else:\n",
    "            mrn2covidcnt[mrn] = covid_count\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean Demo Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:42:42.589051Z",
     "start_time": "2021-12-30T17:42:42.576142Z"
    }
   },
   "outputs": [],
   "source": [
    "#Dictionary to map ethnicity to SOLARStrap friendly coding\n",
    "map_ethnic_dict = { \"NOT HISPANIC OR LATINO OR SPANISH ORIGIN\":'U',\n",
    "\"DECLINED\":'U',\n",
    "\"HISPANIC OR LATINO OR SPANISH ORIGIN\":'H',\n",
    "\"AFRICAN AMERICAN\":'B',\n",
    "\"MULTI-RACIAL\":\"Other\",\n",
    "\"UNKNOWN\":\"U\",\n",
    "\"CAUCASIAN\":\"W\",\n",
    "None:\"U\",\n",
    "\"(null)\":'U',\n",
    "\"ASIAN / PACIFIC ISLANDER\":\"Other\",\n",
    "\"AMERICAN INDIAN / ESKIMO\":\"Other\"\n",
    "                  }\n",
    "\n",
    "map_race_dict = {'WHITE':'W'\n",
    "                    , 'BLACK OR AFRICAN AMERICAN':'B'\n",
    "                    , 'ASIAN':'Other'\n",
    "                    ,'OTHER COMBINATIONS NOT DESCRIBED':'Other'\n",
    "                   ,'AMERICAN INDIAN OR ALASKA NATION':'Other'\n",
    "                   ,'NAT.HAWAIIAN/OTH.PACIFIC ISLAND':'Other'\n",
    "                    ,'SEPHARDIC JEWISH':'Other'\n",
    "                       ,'ASHKENAZI JEWISH':'Other'\n",
    "                    ,'NULL':'NA'\n",
    "                  , '(null)':'NA'\n",
    "                  , 'DECLINED':'NA'\n",
    "              }\n",
    "\n",
    "#Available sex options\n",
    "map_sex_dict = {2:'M',1:'F',3:'NULL',950:'NULL',\n",
    "                \"NULL\":'NULL','(null)':'NULL',\n",
    "                \"M\":\"M\",\"F\":\"F\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:43:12.423440Z",
     "start_time": "2021-12-30T17:43:11.111219Z"
    },
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#Clean race and sex\n",
    "demo_df[\"race_code\"] = demo_df[\"race1\"].map(map_race_dict)\n",
    "demo_df[\"race_code\"].fillna(demo_df[\"race1\"])\n",
    "demo_df[\"sex\"] = demo_df[\"sex_raw\"].map(map_sex_dict)\n",
    "\n",
    "\n",
    "pull_age = lambda x: np.nan if x is None else \\\n",
    "                 (np.nan if x == 0 else date.today().year - x)\n",
    "\n",
    "demo_df[\"age\"] = demo_df[\"birth_year\"].apply(pull_age)\n",
    "\n",
    "demo_df[\"ethnic_code\"] = demo_df[\"ethnicity\"].map(map_ethnic_dict)\n",
    "\n",
    "#GET MISSING empis and add to demo with blank data\n",
    "missing_demo_df = pd.DataFrame( pheno_df[ ~pheno_df[\"ptid\"].\\\n",
    "                                         isin(demo_df[\"mrn\"])][\"ptid\"].values, columns = [\"mrn\"])\n",
    "\n",
    "missing_demo_df.drop_duplicates(inplace=True)\n",
    "for col in demo_df.columns[1:]:\n",
    "    missing_demo_df[col] = None\n",
    "    \n",
    "    \n",
    "#Concat demo with missing demo dataframes\n",
    "demo_all_df = pd.concat([missing_demo_df, demo_df])\n",
    "demo_all_df.reset_index(inplace=True)\n",
    "\n",
    "pull_decade = lambda x: None if x is None\\\n",
    "                 else (None if x =='0' else str(x)[:-1]+\"0\" )\n",
    "demo_all_df[\"birth_decade\"] = demo_all_df[\"birth_year\"].apply(pull_decade)\n",
    "\n",
    "#Fill with \"NA\" for race for SOLAR race/ethnicity Fe algorithm\n",
    "demo_all_df[\"race_code\"].fillna(\"NA\",inplace=True)\n",
    "\n",
    "#Fill with \"NULL\" for other demos if missing\n",
    "fill_with_NULL = [\"age\",\"ethnic_code\",\"sex\"]\n",
    "\n",
    "for demo in fill_with_NULL:\n",
    "    demo_all_df[demo].fillna(\"NULL\",inplace=True)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:43:41.220282Z",
     "start_time": "2021-12-30T17:43:41.206432Z"
    }
   },
   "outputs": [],
   "source": [
    "#Solar implementation on combining race and ethnicity\n",
    "def fe_ethnicity_race(mrn, race, ethnicity):\n",
    "    if race == 'NA':\n",
    "        if ethnicity == 'W':\n",
    "            race = 'White'\n",
    "        elif ethnicity == 'B':\n",
    "            race = 'Black'\n",
    "        elif ethnicity in ('H', 'S', 'O'):\n",
    "            race = 'Hispanic'\n",
    "        elif ethnicity in ('NULL', '', 'U', 'D', '2'):\n",
    "            race = 'Unknown'\n",
    "        else:\n",
    "            race = 'Other'\n",
    "    else:\n",
    "        if race == 'W':\n",
    "            if ethnicity == 'H':\n",
    "                race = 'Hispanic'\n",
    "            else:\n",
    "                race = 'White'\n",
    "        elif race == 'B':\n",
    "            if ethnicity == 'H':\n",
    "                race = 'Hispanic'\n",
    "            else:\n",
    "                race = 'Black'\n",
    "        elif race == 'O':\n",
    "            race = 'Hispanic'\n",
    "        elif race in ('U', 'D'):\n",
    "            if ethnicity == 'H':\n",
    "                race = 'Hispanic'\n",
    "            else:\n",
    "                race = 'Unknown'\n",
    "        else:\n",
    "            race = 'Other'\n",
    "    mrn2race_fe[race][mrn] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T17:44:30.564177Z",
     "start_time": "2021-12-30T17:44:09.795941Z"
    }
   },
   "outputs": [],
   "source": [
    "mrn2race_fe = defaultdict(dict)\n",
    "#Apply fe's alg to rows in df\n",
    "\n",
    "for row in demo_all_df.iterrows():\n",
    "    fe_ethnicity_race(row[1][\"mrn\"], row[1][\"race_code\"],row[1][\"ethnic_code\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T18:09:10.379870Z",
     "start_time": "2021-12-30T18:09:10.371322Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "202798"
      ]
     },
     "execution_count": 344,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(demo_all_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-30T18:08:40.784824Z",
     "start_time": "2021-12-30T18:08:39.916215Z"
    }
   },
   "outputs": [],
   "source": [
    "#Output Demo file\n",
    "file_name = filename_demo\n",
    "\n",
    "demo_all_df[[\"mrn\",\"sex\",\"birth_decade\",\"race_code\",\"ethnic_code\",\"age\"]].\\\n",
    "to_csv(file_name,index=False, sep='\\t', header=[\"ptid\",\"sex\",\"birth_decade\",\"race_code\",\n",
    "                                           \"ethnic_code\",\"age\"])\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheno: Cummulative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-13T21:50:00.846913Z",
     "start_time": "2022-01-13T21:50:00.835942Z"
    }
   },
   "outputs": [],
   "source": [
    "def output_df_function(df, filename_or_pheno, base_path):\n",
    "#Output covid file\n",
    "    file_name = base_path + filename_or_pheno + \".txt\"\n",
    "    \n",
    "    if not os.path.isfile(file_name):\n",
    "        df.to_csv(file_name,index=False, sep='\\t',\\\n",
    "                  mode='a',columns=[\"ptid\",\"pheno\",\"value\"] )\n",
    "   \n",
    "    else: # else it exists so append without writing the header\n",
    "        df.to_csv(file_name,index=False, sep='\\t',\\\n",
    "                   mode='a', header=False, columns=[\"ptid\",\"pheno\",\"value\"] )\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-11T21:52:43.467276Z",
     "start_time": "2022-02-11T21:52:43.446939Z"
    },
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#Get list of pheno datecut dataframes to concat\n",
    "\n",
    "def pull_datecut_dataframe(df, allow_class_switching=False):\n",
    "    \n",
    "    #Get earliest and latest test dates\n",
    "    date = min( df[\"order_date\"] )\n",
    "    last_test_date = max( df[\"order_date\"] )\n",
    "    \n",
    "    #All phenos \n",
    "    pheno_datecut_list = list()\n",
    "    \n",
    "    #Dates of phenos\n",
    "    datecut_date_list = list()\n",
    "    \n",
    "    #Number of families with those phenos\n",
    "    datecut_fam_num_list = list()\n",
    "    \n",
    "    naming_pheno = \"covid_weekly_\"\n",
    "\n",
    "    while date < last_test_date:\n",
    "        date += timedelta(days= 7)\n",
    "        \n",
    "        month = date.month\n",
    "        year = date.year\n",
    "\n",
    "        #get all first positive (for cases)\n",
    "        case_df = df[ ( df[\"value\"] == 1 ) &\\\n",
    "                        (df['order_date'] <= date)] \n",
    "        #get all controls irrespective of date cutoff - assume always a case if tested\n",
    "        control_df = df_control.copy() \n",
    "        \n",
    "        if allow_class_switching:\n",
    "            naming_pheno = \"covid_cswitch_weekly_\"\n",
    "            \n",
    "            #These are patients that will be infected but jhavente been yet\n",
    "            future_cases_df = df[ ( df[\"value\"] == 1 ) &\\\n",
    "                (df['order_date'] > date)] \n",
    "            \n",
    "            control_df = pd.concat([future_cases_df, control_df])        \n",
    "            \n",
    "            #Turn future cases into controls\n",
    "            control_df[\"value\"] = 0\n",
    "\n",
    "        df_datecut = pd.concat([case_df, control_df])\n",
    "        \n",
    "        pheno = naming_pheno + date.strftime('%m_%d_%Y')\n",
    "        df_datecut[\"pheno\"] = pheno\n",
    "        \n",
    "        #Get the number of SOLAR families that qualify with that date cut (ie >=1 case and >1 case+control)\n",
    "        fam_num = pull_fam_count(pheno, df_datecut, print_res = False)\n",
    "        \n",
    "        #Write the sample sizes\n",
    "        write_famcount( fam_num, naming_pheno + date.strftime('%m_%d_%Y'))\n",
    "        \n",
    "        datecut_fam_num_list.append(fam_num)\n",
    "        datecut_date_list.append( date )\n",
    "        \n",
    "        #Set minimum number of solar families to be 45\n",
    "        if fam_num >= 45:\n",
    "            pheno_datecut_list.append(df_datecut) \n",
    "            filename = naming_pheno +\\\n",
    "                        str( month ) +\"_\" + str( year)\n",
    "            \n",
    "            base_path = '/home/kel2158/Projects/COVIDHerit/data/12_29_21/cummulative/'\n",
    "            \n",
    "            if allow_class_switching:\n",
    "                base_path = '/home/kel2158/Projects/COVIDHerit/data/12_29_21/cummulative_classswitch/'\n",
    "            \n",
    "            output_df_function(df_datecut, filename, base_path)\n",
    "            \n",
    "#     Concat \n",
    "    pheno_datecut_df = pd.concat(pheno_datecut_list)\n",
    "    \n",
    "    return( pheno_datecut_df )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-28T18:06:04.480819Z",
     "start_time": "2022-01-28T18:05:48.549960Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Create control dataframe\n",
    "df_control = pd.DataFrame(controls,columns=[\"ptid\"])\n",
    "df_control[\"value\"] = 0\n",
    "df_control[\"order_time\"] = np.nan\n",
    "df_control[\"famid\"] = df_control[\"ptid\"].map(mrn2famid)\n",
    "\n",
    "#Pull dataframe\n",
    "pheno_datecut_df = pull_datecut_dataframe(first_pcr_nohospreq_df)\n",
    "\n",
    "#Pull dataframe\n",
    "pheno_classswitch_datecut_df = pull_datecut_dataframe(first_pcr_nohospreq_df, allow_class_switching=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-11T21:57:00.045427Z",
     "start_time": "2022-02-11T21:56:57.940886Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of weeks with over 40 families when a control cannot switch to a case: 84 \n",
      "First week with over 40 families: covid_weekly_03_24_2020 \n"
     ]
    }
   ],
   "source": [
    "print( f\"Number of weeks with over 40 families when a control cannot switch to a case: {len( set( pheno_datecut_df['pheno']  ) )} \" )\n",
    "      \n",
    "print( f\"First week with over 40 families: { pheno_datecut_df['pheno'].values[0] } \" )                                                         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-11T21:56:24.558892Z",
     "start_time": "2022-02-11T21:56:22.380338Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of weeks with over 40 families when a control can switch to a case: 84 \n",
      "First week with over 40 families: covid_cswitch_weekly_03_24_2020 \n"
     ]
    }
   ],
   "source": [
    "\n",
    "print( f\"Number of weeks with over 40 families when a control can switch to a case: {len( set( pheno_classswitch_datecut_df['pheno']  ) )} \" )\n",
    "      \n",
    "print( f\"First week with over 40 families: { pheno_classswitch_datecut_df['pheno'].values[0] } \" )                                                         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheno: Days Hospitalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 354,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:06:24.144572Z",
     "start_time": "2021-12-31T16:06:24.136075Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_days(timedelta_date):\n",
    "\n",
    "    days = timedelta_date.days\n",
    "    partial_day = timedelta_date.seconds/(24*60*60)\n",
    "    \n",
    "    return(days+partial_day)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 355,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:08:03.695775Z",
     "start_time": "2021-12-31T16:06:59.640432Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 11979/11979 [01:03<00:00, 187.92it/s]\n"
     ]
    }
   ],
   "source": [
    "#We want to limit first hosp dataframe so every patient has just one initial start encounter for getting covid test\n",
    "#Select on earliest admission time\n",
    "#Why not just use defined_visits?  Some are wonky - ex pat_enc_csn_id = 147040327\n",
    "grouped_firstadmin = first_admsn_df[[\"pat_enc\",\"mrn\",\"admsn_time\",\"disch_time\"]].drop_duplicates().groupby([\"mrn\"])\n",
    "\n",
    "mrn2first_hospenc = dict()\n",
    "mrn2first_hospenc_admsn = dict()\n",
    "mrn2first_hospenc_disch = dict()\n",
    "\n",
    "for group in tqdm( grouped_firstadmin):\n",
    "\n",
    "    mrn = group[0]\n",
    "    holder_df = group[1]\n",
    "    \n",
    "    #First filter by getting the earliest admission within 96 hours of first PCR test\n",
    "    holder_df2 = holder_df[ holder_df[\"admsn_time\"] == holder_df[\"admsn_time\"].min() ] \n",
    "\n",
    "    #Then filter by getting the latest discharge time associated with that intake time\n",
    "    #ie may transfer between centers etc\n",
    "    holder_df3 = holder_df2[ holder_df2[\"disch_time\"] == holder_df2[\"disch_time\"].max() ] \n",
    "\n",
    "    #Sometimes identical patient encounters in terms of admission and discharge time (Different department location) \n",
    "    #Don't care about this so just take first encounter\n",
    "    patient_encounter = holder_df3[\"pat_enc\"].values[0]\n",
    "    patient_admsn = holder_df3[\"admsn_time\"].values[0]\n",
    "    patient_disch = holder_df3[\"disch_time\"].values[0]\n",
    "    \n",
    "\n",
    "    #Sometimes the admission time and discharge time is slightly different for a single encounter.\n",
    "    mrn2first_hospenc[mrn] = patient_encounter\n",
    "    mrn2first_hospenc_admsn[mrn] = patient_admsn\n",
    "    mrn2first_hospenc_disch[mrn] = patient_disch\n",
    "\n",
    "    \n",
    "    #Making sure they are all in fact identical except for patient encounter number\n",
    "    if len(holder_df3[[\"mrn\",\"admsn_time\",\"disch_time\"]].drop_duplicates()) >1:\n",
    "        print(holder_df3)\n",
    "        break\n",
    "\n",
    "\n",
    "#Limit first hosp dataframe so every patient has just one initial start encounter\n",
    "first_admsn_limit_df = first_admsn_df[ ( first_admsn_df[\"mrn\"].map(mrn2first_hospenc) == first_admsn_df[\"pat_enc\"]) \n",
    "                                     & ( first_admsn_df[\"mrn\"].map(mrn2first_hospenc_admsn) == first_admsn_df[\"admsn_time\"])\n",
    "                                     & ( first_admsn_df[\"mrn\"].map(mrn2first_hospenc_disch) == first_admsn_df[\"disch_time\"]) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 356,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:08:39.505599Z",
     "start_time": "2021-12-31T16:08:39.426313Z"
    }
   },
   "outputs": [],
   "source": [
    "#Drop duplicates\n",
    "first_admsn_limit_df.drop(columns=\"visit_type\",inplace=True,axis=1)\n",
    "first_admsn_limit_df.drop_duplicates(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 357,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:09:14.451530Z",
     "start_time": "2021-12-31T16:09:14.437255Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 357,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Each mrn should have one first admission tied to a positive test\n",
    "len(first_admsn_limit_df) == len(set(first_admsn_limit_df[\"mrn\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-29T19:36:18.617334Z",
     "start_time": "2021-09-29T19:36:18.507230Z"
    }
   },
   "outputs": [],
   "source": [
    "#Limit admissions to those that are discharged\n",
    "admsn_df = discharged_admissions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:12:19.389249Z",
     "start_time": "2021-12-31T16:12:19.372793Z"
    }
   },
   "outputs": [],
   "source": [
    "#Identical patient encounter IDs with admission, transfer etc for visit type.  To drop visit_type and dedup\n",
    "\n",
    "first_admsn_limit_df.drop_duplicates(inplace=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:13:31.324817Z",
     "start_time": "2021-12-31T16:13:29.064593Z"
    }
   },
   "outputs": [],
   "source": [
    "#Get time in hospital\n",
    "admsn_df[\"hosp_days\"] = admsn_df[\"hosp_disch\"] - admsn_df[\"hosp_admsn\"]\n",
    "admsn_df[\"hosp_days\"] = admsn_df[\"hosp_days\"].apply(get_days)\n",
    "\n",
    "#Sort values by admission time and then mrn\n",
    "admsn_df.sort_values([\"hosp_admsn\"],ascending= True,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:14:26.602000Z",
     "start_time": "2021-12-31T16:14:06.625492Z"
    }
   },
   "outputs": [],
   "source": [
    "mrn2admsissions = defaultdict(dict)\n",
    "\n",
    "for mrn in first_admsn_limit_df[\"mrn\"]:\n",
    "    #Only include those admissions that last longer than a day, otherwise don't care\n",
    "    df = admsn_df[ (admsn_df[\"mrn\"] == mrn) & (admsn_df[\"hosp_days\"] >= 1) ]\n",
    "    adm_lst = df[\"hosp_admsn\"].to_list()\n",
    "    disch_lst = df[\"hosp_disch\"].to_list()\n",
    "    \n",
    "    mrn2admsissions[mrn][\"admsn\"] = adm_lst\n",
    "    mrn2admsissions[mrn][\"disch\"] = disch_lst\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:16:22.829956Z",
     "start_time": "2021-12-31T16:16:20.535101Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "11979it [00:02, 5278.53it/s]\n"
     ]
    }
   ],
   "source": [
    "#Loop through and see if re-admitted for covid\n",
    "#ie if never admitted, admitted within 12 days\n",
    "#if admitted, admitted within 5 days of previous discharge\n",
    "\n",
    "mrn2final_disch = dict()\n",
    "mrn2final_admission = dict()\n",
    "mrn2hosp_days = dict()\n",
    "mrn2already_hosp = dict()\n",
    "mrn_openadmission = list()\n",
    "\n",
    "for row in tqdm( first_admsn_limit_df.iterrows() ):\n",
    "    mrn = row[1][\"mrn\"]\n",
    "    hosp_days = float( row[1][\"hrs_from_admission\"] )\n",
    "    initial_admsn_date = row[1][\"admsn_time\"]\n",
    "    initial_disch_date = row[1][\"disch_time\"]\n",
    "    \n",
    "    \n",
    "    #Assume only 1 admission/discharge related to covid start\n",
    "    final_admsn = initial_admsn_date\n",
    "    final_disch = initial_disch_date\n",
    "    \n",
    "    #Note, these are limited to admissions lasting more than 24 hours\n",
    "    further_adm = mrn2admsissions[mrn][\"admsn\"]\n",
    "    further_dis = mrn2admsissions[mrn][\"disch\"]\n",
    "    \n",
    "    #Reset hosp_days to 0,admitted an inconsequential amount of time (under a day) or outpatient\n",
    "    if hosp_days < 1:\n",
    "        hosp_days = 0\n",
    "    \n",
    "    num_admin = len(further_adm)\n",
    "    \n",
    "    \n",
    "    #If patient has been in there less than a day\n",
    "    #see if there is another admission within 12 days of initial discharge that lasted more than a day \n",
    "    if num_admin > 0:\n",
    "        \n",
    "        for admission in zip( further_adm,further_dis):\n",
    "            nxt_admsn = admission[0]\n",
    "            nxt_disch = admission[1]\n",
    "            \n",
    "            #Admission is before initial covid admission\n",
    "            if nxt_admsn < initial_admsn_date:\n",
    "                \n",
    "                #Need to make sure patient wasn't already in hospital when they got the covid test\n",
    "                if nxt_disch >= initial_admsn_date:\n",
    "                    \n",
    "                    #Problem, looks like patient was already in hospital prior to covid test and may have gotten it \n",
    "                    #There\n",
    "                    mrn2already_hosp[mrn]  = [nxt_admsn,nxt_disch]\n",
    " \n",
    "            #Re-admitted within 12 days of initial positive covid pcr test which had been outpatient/ER\n",
    "            #ie initial hospital stay w/ covid was for less than a day\n",
    "            elif hosp_days == 0:\n",
    "\n",
    "                #Patient declines within 12 days of first positive and admitted (ie there at least a day)\n",
    "                if get_days(nxt_admsn - final_disch) <= 12:\n",
    "\n",
    "                    #Add additional days to total days in hospital\n",
    "                    hosp_days += get_days(nxt_disch - nxt_admsn)\n",
    "       \n",
    "\n",
    "                    #Update final discharge date\n",
    "                    final_disch = nxt_disch\n",
    "                    final_admsn = nxt_admsn\n",
    "                \n",
    "                \n",
    "            #There are sometimes overlapping admissions, see if current admission occurs before final discharge\n",
    "            elif nxt_admsn < final_disch:\n",
    "                \n",
    "                #Check if subsequent discharge is later\n",
    "                if nxt_disch > final_disch:\n",
    "                    \n",
    "                    #If so, need to add differences between discharges and update final discharge\n",
    "                    hosp_days += get_days(nxt_disch - final_disch)\n",
    "                    \n",
    "                    #Assign this later discharge as final\n",
    "                    final_disch = nxt_disch\n",
    "                    \n",
    "                    \n",
    "            #Else patient was initially admitted, see if readmitted b/c of covid (ie within 5 days of previous)\n",
    "            #This can be negative if overlapping admissions/discharges\n",
    "            elif get_days(nxt_admsn - final_disch) <= 5:\n",
    "                hosp_days += get_days(nxt_disch - nxt_admsn)\n",
    "                final_disch = nxt_disch\n",
    "                final_admsn = nxt_admsn\n",
    "                \n",
    "    mrn2final_disch[mrn] = final_disch\n",
    "    mrn2final_admission[mrn] = final_admsn\n",
    "    mrn2hosp_days[mrn] = hosp_days\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 374,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:22:12.379667Z",
     "start_time": "2021-12-31T16:22:10.971984Z"
    }
   },
   "outputs": [],
   "source": [
    "first_admsn_limit_df[\"hosp_days\"] = first_admsn_limit_df[\"mrn\"].map(mrn2hosp_days)\n",
    "first_admsn_limit_df[\"final_disch\"] = first_admsn_limit_df[\"mrn\"].map(mrn2final_disch)\n",
    "\n",
    "#These are people that were hospitalized during time of covid test\n",
    "first_admsn_limit_df[\"already_hosp\"] = first_admsn_limit_df[\"mrn\"].map(mrn2already_hosp)\n",
    "\n",
    "#Map people that have died\n",
    "mrn2death = dict(zip( deaths_df[\"mrn\"],deaths_df[\"death_date\"]))\n",
    "first_admsn_limit_df[\"death_date\"] = first_admsn_limit_df[\"mrn\"].map(mrn2death)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:28:32.323892Z",
     "start_time": "2021-12-31T16:28:32.231733Z"
    }
   },
   "outputs": [],
   "source": [
    "#limit to patients that were not already in hospital prior to first positive test (ie over 96 hrs)\n",
    "#And no open admissions\n",
    "#And their final discharge is at least 12 days before data pull \n",
    "#And patient hasn't died within 5 days of discharge\n",
    "\n",
    "data_pull = first_admsn_limit_df[\"final_disch\"].max()\n",
    "\n",
    "hosp_patient_df = first_admsn_limit_df[ ( first_admsn_limit_df[\"already_hosp\"].isnull()) &\n",
    "                    ( ~ first_admsn_limit_df[\"mrn\"].isin( open_admission[\"mrn\"] )) & \n",
    "                    ( (first_admsn_limit_df[\"final_disch\"] + timedelta(12)) < data_pull ) &\n",
    "                    ( ( first_admsn_limit_df[\"death_date\"].isnull())  | \\\n",
    "                        ( first_admsn_limit_df[\"death_date\"] - timedelta(5)  >  first_admsn_limit_df[\"final_disch\"]))  \n",
    "                                    \n",
    "                                      ] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosp_patient_df[\"famid\"] = hosp_patient_df[\"mrn\"].map(mrn2famid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T17:40:38.399549Z",
     "start_time": "2021-12-31T17:40:38.387384Z"
    }
   },
   "outputs": [],
   "source": [
    "pheno_days_hosp = hosp_patient_df[[ \"mrn\",\"hosp_days\",\"famid\" ]]\n",
    "pheno_days_hosp.rename(columns={\"mrn\":\"ptid\",\"hosp_days\":\"value\"},inplace=True)\n",
    "pheno_days_hosp[\"pheno\"] = \"days_hosp\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 432,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T17:56:04.951862Z",
     "start_time": "2021-12-31T17:56:04.853276Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 115 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"days_hosp\", pheno_days_hosp),\\\n",
    "               \"days_hosp\")\n",
    "\n",
    "#Output Outcomes\n",
    "file_name = output_folder + \"days_hosp.txt\"\n",
    "\n",
    "pheno_days_hosp[cols_to_keep].to_csv(file_name,index=False, sep='\\t')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-09-29T19:36:40.543995Z",
     "start_time": "2021-09-29T19:36:40.481123Z"
    }
   },
   "source": [
    "## Pheno: Intubation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:38:41.532259Z",
     "start_time": "2021-12-31T16:38:41.493387Z"
    }
   },
   "outputs": [],
   "source": [
    "#Combine admissions data with intubation data\n",
    "#Dont care if patient has died or discharged or whatever\n",
    "intub_hospcombined_df = intub_df.merge( first_admsn_limit_df,\n",
    "                                  on=\"mrn\")\n",
    "\n",
    "#Make datetime (not time excluded)\n",
    "intub_hospcombined_df[\"intubation_date\"] = pd.to_datetime(intub_hospcombined_df[\"intubation_date\"])\n",
    "\n",
    "#Patients intubated after covid admission start (add a day since no time) with test, \n",
    "#but before discharge\n",
    "intub_covidstay_df = intub_hospcombined_df[ \n",
    "    ( intub_hospcombined_df[\"intubation_date\"] + timedelta(1) > intub_hospcombined_df[\"admsn_time\"] ) &\n",
    "( intub_hospcombined_df[\"intubation_date\"] < intub_hospcombined_df[\"final_disch\"] )  ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:39:24.962635Z",
     "start_time": "2021-12-31T16:39:24.838359Z"
    }
   },
   "outputs": [],
   "source": [
    "#Cases are patients that have been intubated during a covid stay\n",
    "mrn2intub_covidstay = dict(zip( intub_covidstay_df.mrn.unique(), \n",
    "                               len(intub_covidstay_df.mrn.unique())*[1]) )\n",
    "#Controls are all patients with positive pcr test that don't have a dni/dnr\n",
    "\n",
    "intub_pheno_df = pheno_df[pheno_df[\"value\"] == 1][[\"ptid\",\"famid\"]]\n",
    "#Being very conservative, irrespective of date/covid date, exclude if there is a dnr/dni order\n",
    "intub_pheno_df[\"dni_dnr_date\"] = intub_pheno_df[\"ptid\"].map(mrn2dni_dnr)\n",
    "intub_pheno_df[\"value\"] = intub_pheno_df[\"ptid\"].map(mrn2intub_covidstay)\n",
    "intub_pheno_df[\"value\"].fillna(0,inplace=True)\n",
    "intub_pheno_df[\"pheno\"] = \"intubation\"\n",
    "\n",
    "#If hospitalized include final discharge\n",
    "intub_pheno_df[\"final_disch_ifhosp\"] =intub_pheno_df[\"ptid\"].map(mrn2final_disch)\n",
    "\n",
    "#Patients with DNI/DNR specified before discharge from covid stay\n",
    "#and were intubated during covid stay\n",
    "dnr_dni_duringcovidstay_df = intub_pheno_df[ ( intub_pheno_df[\"dni_dnr_date\"].notnull() )& \n",
    "(intub_pheno_df[\"final_disch_ifhosp\"] > pd.to_datetime( intub_pheno_df[\"dni_dnr_date\"])) \n",
    "              ]\n",
    "#Final dataframe, exclude patients with DNRs\n",
    "intub_pheno_df = intub_pheno_df[~ (\n",
    "    ( intub_pheno_df[\"ptid\"].isin( dnr_dni_duringcovidstay_df[\"ptid\"].unique() ) )\n",
    "    & (intub_pheno_df[\"ptid\"] == 0 ) \n",
    "                                    )]\n",
    "\n",
    "intub_pheno_df.drop_duplicates(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-06T15:24:42.960623Z",
     "start_time": "2022-01-06T15:24:42.937599Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 26 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"intubation\",intub_pheno_df),\\\n",
    "               \"intubation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheno Death"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:40:07.883275Z",
     "start_time": "2021-12-31T16:40:06.687817Z"
    }
   },
   "outputs": [],
   "source": [
    "#Cases are patients with positive covid test that die\n",
    "#Controls are those that don;t die\n",
    "\n",
    "death_pheno_df = pheno_df[pheno_df[\"value\"] == 1][[\"ptid\",\"famid\"]]\n",
    "#Being very conservative, irrespective of date/covid date, exclude if there is a dnr/dni order\n",
    "death_pheno_df[\"dni_dnr_date\"] = death_pheno_df[\"ptid\"].map(mrn2dni_dnr)\n",
    "death_pheno_df[\"death_datetime\"] = death_pheno_df[\"ptid\"].map(mrn2death)\n",
    "\n",
    "#If hospitalized include final discharge\n",
    "death_pheno_df[\"final_disch_ifhosp\"] =death_pheno_df[\"ptid\"].map(mrn2final_disch)\n",
    "death_pheno_df[\"pheno\"] = \"death\"\n",
    "\n",
    "\n",
    "#Get patients who died within 30 days of covid hospitalization stay\n",
    "death_covidhosp_df = death_pheno_df [ (death_pheno_df[\"death_datetime\"].notnull()) & \n",
    "(death_pheno_df[\"death_datetime\"]  < death_pheno_df[\"final_disch_ifhosp\"]+ timedelta(30) ) ]\n",
    "\n",
    "mrn2death_duringcovhosp = dict(zip( death_covidhosp_df[\"ptid\"],len( death_covidhosp_df[\"ptid\"])*[1] ))\n",
    "\n",
    "\n",
    "#For now, don't care about DNI/DNR order\n",
    "death_pheno_df[\"value\"] = death_pheno_df[\"ptid\"].map(mrn2death_duringcovhosp)\n",
    "death_pheno_df[\"value\"].fillna(0,inplace=True)\n",
    "\n",
    "\n",
    "death_pheno_df.drop_duplicates(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-06T15:24:56.637373Z",
     "start_time": "2022-01-06T15:24:56.617187Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 0 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"death\",death_pheno_df),\\\n",
    "               \"death\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pheno: Hospitalization Status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:40:49.990677Z",
     "start_time": "2021-12-31T16:40:49.951419Z"
    }
   },
   "outputs": [],
   "source": [
    "#Hospital status where 0/1 for hospitalized or not\n",
    "mrn2hosp_status = dict(zip( hosp_patient_df[ hosp_patient_df[\"hosp_days\"] == 0][\"mrn\"] ,\n",
    "         len(hosp_patient_df[ hosp_patient_df[\"hosp_days\"] < 1][\"mrn\"])*[0]))\n",
    "\n",
    "mrn2hosp_status.update( dict(zip( hosp_patient_df[ hosp_patient_df[\"hosp_days\"] > 0][\"mrn\"] ,\n",
    "         len(hosp_patient_df[ hosp_patient_df[\"hosp_days\"] >= 1][\"mrn\"])*[1])))\n",
    "\n",
    "\n",
    "hosp_patient_df[\"hosp_stat\"] = hosp_patient_df[\"mrn\"].map(mrn2hosp_status)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:41:31.846264Z",
     "start_time": "2021-12-31T16:41:31.829043Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6612027972027972"
      ]
     },
     "execution_count": 385,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Positive covid patients we have hosp info for\n",
    "len(hosp_patient_df.mrn.unique())/len( pheno_df[ pheno_df[\"value\"] == 1 ][\"ptid\"].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T16:42:14.613351Z",
     "start_time": "2021-12-31T16:42:14.593196Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6701538461538461"
      ]
     },
     "execution_count": 386,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Large falloff when we limit to those patients that get a covid test within\n",
    "#96 hrs of admission (no limit in length of time)\n",
    "\n",
    "#Additional falloff with patients excluded inculde (filtered out eaerlier in script): \n",
    "#1. have open admission, \n",
    "#2. are discharged within 12 days of data pull\n",
    "#3. have died\n",
    "\n",
    "len(first_admsn_df.mrn.unique())/len( pheno_df[ pheno_df[\"value\"] == 1 ][\"ptid\"].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 407,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T17:18:28.218835Z",
     "start_time": "2021-12-31T17:18:28.203825Z"
    }
   },
   "outputs": [],
   "source": [
    "pheno_hospstat_df = hosp_patient_df[hosp_patient_df[\"hosp_stat\"].notnull()][[\"mrn\",\"hosp_stat\",\"famid\"]]\n",
    "\n",
    "#Format\n",
    "pheno_hospstat_df[\"pheno\"] = \"hosp_stat\"\n",
    "pheno_hospstat_df.rename(columns={\"mrn\": \"ptid\", \"hosp_stat\":\"value\" },inplace=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 413,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T17:24:23.295954Z",
     "start_time": "2021-12-31T17:24:23.274911Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 115 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"hosp_stat\", pheno_hospstat_df),\\\n",
    "               \"hosp_stat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T18:02:58.833204Z",
     "start_time": "2021-12-31T18:02:58.747130Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 115 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "write_famcount( pull_fam_count(\"hosp_stat\", pheno_hospstat_df),\\\n",
    "               \"hosp_stat\")\n",
    "\n",
    "#Output hosp stat outcome\n",
    "file_name = output_folder + \"hosp_stat.txt\"\n",
    "\n",
    "#pheno_binary_df.to_csv(file,index=False, sep='\\t')\n",
    "pheno_hospstat_df[cols_to_keep].to_csv(file_name,index=False, sep='\\t')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 447,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-03T18:49:03.751503Z",
     "start_time": "2022-01-03T18:49:03.725502Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 158 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "pheno_hospstat_df[\"pheno\"] = \"hosp_stat_nopro\"\n",
    "\n",
    "write_famcount( pull_fam_count(\"hosp_stat_nopro\", pheno_hospstat_df,type_pheno=\"noproband\" ),\\\n",
    "               \"hosp_stat_nopro\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-31T18:01:45.768412Z",
     "start_time": "2021-12-31T18:01:45.716594Z"
    }
   },
   "outputs": [],
   "source": [
    "#Output with no proband\n",
    "file_name = output_folder + \"hosp_stat_noproband.txt\"\n",
    "\n",
    "pheno_hospstat_df[cols_to_keep].to_csv(file_name,index=False, sep='\\t')\n",
    "\n",
    "#Go back to normal\n",
    "pheno_hospstat_df[\"pheno\"] = \"hosp_stat\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "covid_herit",
   "language": "python",
   "name": "covid_herit"
  },
  "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.10"
  },
  "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": false,
   "toc_position": {
    "height": "356px",
    "left": "59px",
    "top": "494px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}