Tables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c41400c3",
   "metadata": {},
   "source": [
    "\n",
    "--* Table1: 3-way ANOVA\n",
    "* Table2: 2-way anova, alberta data\n",
    "* Table S1: Genome list, accessions, habitat, clade A/B, other metadata.\n",
    "* Table S3: Presence/Absence of pan-genome and targeted genes by genome\n",
    "* Table S5: Pagel data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2c1201a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "import scikit_posthocs as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2db4dc62",
   "metadata": {},
   "outputs": [],
   "source": [
    "long = pd.read_table('genomes_vs_features_longform.csv', sep=',')\n",
    "meta = pd.read_table('genome_labels.csv', sep=',')\n",
    "\n",
    "def rename_clade(string):\n",
    "    if type(string) == str:\n",
    "        return string.replace('Clade', 'Type')\n",
    "    return string\n",
    "\n",
    "long['Type'] = long['Clade'].map(rename_clade)\n",
    "loc_abbr_map = {'United Kingdom': 'UK', \n",
    "                'Canada/Alberta': 'Ab'}\n",
    "long['Sampling Location'] = long.apply(lambda row: ' '.join([row['Habitat'], loc_abbr_map[row['Country/Province']]]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "36d21a34",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = long[['Isolate', 'Country/Province', 'Habitat', 'Type', 'source', 'Presence']]\\\n",
    "            .groupby(['Isolate', 'Country/Province', 'Habitat', 'Type', 'source',])\\\n",
    "            .sum().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7debcb99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Isolate</th>\n",
       "      <th>Country/Province</th>\n",
       "      <th>Habitat</th>\n",
       "      <th>Type</th>\n",
       "      <th>source</th>\n",
       "      <th>Presence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>AMR</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Genomic Island</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Metal</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Novel Plasmid</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Phage</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8899</th>\n",
       "      <td>SRR14026555</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clinical</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Metal</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8900</th>\n",
       "      <td>SRR14026555</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clinical</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Novel Plasmid</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8901</th>\n",
       "      <td>SRR14026555</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clinical</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Phage</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8902</th>\n",
       "      <td>SRR14026555</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clinical</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8903</th>\n",
       "      <td>SRR14026555</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clinical</td>\n",
       "      <td>Type A</td>\n",
       "      <td>VF</td>\n",
       "      <td>35.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8904 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          Isolate Country/Province       Habitat    Type          source  \\\n",
       "0      ERR1007500   United Kingdom  Agricultural  Type A             AMR   \n",
       "1      ERR1007500   United Kingdom  Agricultural  Type A  Genomic Island   \n",
       "2      ERR1007500   United Kingdom  Agricultural  Type A           Metal   \n",
       "3      ERR1007500   United Kingdom  Agricultural  Type A   Novel Plasmid   \n",
       "4      ERR1007500   United Kingdom  Agricultural  Type A           Phage   \n",
       "...           ...              ...           ...     ...             ...   \n",
       "8899  SRR14026555   Canada/Alberta      Clinical  Type A           Metal   \n",
       "8900  SRR14026555   Canada/Alberta      Clinical  Type A   Novel Plasmid   \n",
       "8901  SRR14026555   Canada/Alberta      Clinical  Type A           Phage   \n",
       "8902  SRR14026555   Canada/Alberta      Clinical  Type A         Plasmid   \n",
       "8903  SRR14026555   Canada/Alberta      Clinical  Type A              VF   \n",
       "\n",
       "      Presence  \n",
       "0         18.0  \n",
       "1          4.0  \n",
       "2          5.0  \n",
       "3          0.0  \n",
       "4          3.0  \n",
       "...        ...  \n",
       "8899       6.0  \n",
       "8900       0.0  \n",
       "8901       4.0  \n",
       "8902       5.0  \n",
       "8903      35.0  \n",
       "\n",
       "[8904 rows x 6 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7332087d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Isolate</th>\n",
       "      <th>Feature</th>\n",
       "      <th>Presence</th>\n",
       "      <th>Habitat</th>\n",
       "      <th>Country/Province</th>\n",
       "      <th>Clade</th>\n",
       "      <th>source</th>\n",
       "      <th>Type</th>\n",
       "      <th>Sampling Location</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SRR14026455</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Wastewater Agr.</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade B</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Type B</td>\n",
       "      <td>Wastewater Agr. Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SRR14026462</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SRR14026496</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SRR14026467</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SRR14026460</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade B</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Type B</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Isolate Feature  Presence          Habitat Country/Province    Clade  \\\n",
       "0  SRR14026455   AA083       0.0  Wastewater Agr.   Canada/Alberta  Clade B   \n",
       "1  SRR14026462   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "2  SRR14026496   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "3  SRR14026467   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "4  SRR14026460   AA083       0.0    Natural Water   Canada/Alberta  Clade B   \n",
       "\n",
       "    source    Type   Sampling Location  \n",
       "0  Plasmid  Type B  Wastewater Agr. Ab  \n",
       "1  Plasmid  Type A    Natural Water Ab  \n",
       "2  Plasmid  Type A    Natural Water Ab  \n",
       "3  Plasmid  Type A    Natural Water Ab  \n",
       "4  Plasmid  Type B    Natural Water Ab  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "long.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a19a13b",
   "metadata": {},
   "source": [
    "## Table X: Unassigned. Feature presences/percentages by sampling location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "64e2c8ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "for feat in long['source'].unique():\n",
    "    data = long[long['source']==feat]\n",
    "    sort_order = ['Presence', 'Percentage']\n",
    "    frames = []\n",
    "    \n",
    "    total = long[long['source']==feat].groupby('Feature').sum()\n",
    "    total = total.join(total.div(1273).rename({'Presence': \"Percentage\"}, axis=1))\n",
    "    total.sort_values(by='Presence', ascending=False, inplace=True)\n",
    "    \n",
    "    for col in ['Type', 'Sampling Location']:\n",
    "        p = data.groupby([col, 'Feature']).agg({'Presence': \"sum\"})\n",
    "        c = data.groupby([col]).agg({'Isolate': \"nunique\"})\n",
    "        c.rename({'Isolate': \"Presence\"}, axis=1, inplace=True)\n",
    "        perc = p.div(c, level=col).reset_index().pivot(index='Feature', columns=col, values='Presence')\n",
    "        presence = p.reset_index().pivot(index='Feature', columns=col, values='Presence')\n",
    "\n",
    "        perc.columns.name = ''\n",
    "        presence.columns.name =''\n",
    "\n",
    "        presence.rename({column: '{} Number Present'.format(column) for column in presence.columns}, axis=1, inplace=True)\n",
    "        perc.rename({column: '{} Percent Present'.format(column) for column in perc.columns}, axis=1, inplace=True)\n",
    "\n",
    "        vals = sorted(data[col].dropna().unique())\n",
    "        for val in vals:\n",
    "            sort_order.append('{} Number Present'.format(val))\n",
    "            sort_order.append('{} Percent Present'.format(val))\n",
    "\n",
    "        frames.append(presence.join(perc))\n",
    "    total.join(frames)[sort_order].to_csv('tables/feature_presence_percentages/{}_summary.csv'.format(feat), sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "337bd041",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Isolate</th>\n",
       "      <th>Country/Province</th>\n",
       "      <th>Habitat</th>\n",
       "      <th>Type</th>\n",
       "      <th>source</th>\n",
       "      <th>Presence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>AMR</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Genomic Island</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Metal</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Novel Plasmid</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Phage</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Isolate Country/Province       Habitat    Type          source  Presence\n",
       "0  ERR1007500   United Kingdom  Agricultural  Type A             AMR      18.0\n",
       "1  ERR1007500   United Kingdom  Agricultural  Type A  Genomic Island       4.0\n",
       "2  ERR1007500   United Kingdom  Agricultural  Type A           Metal       5.0\n",
       "3  ERR1007500   United Kingdom  Agricultural  Type A   Novel Plasmid       0.0\n",
       "4  ERR1007500   United Kingdom  Agricultural  Type A           Phage       3.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc068eab",
   "metadata": {},
   "source": [
    "## Table 1: 3-way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4a61539c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alex/anaconda3/envs/biopython/lib/python3.9/site-packages/pandas/core/frame.py:5039: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  return super().rename(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Isolate</th>\n",
       "      <th>Geography</th>\n",
       "      <th>Habitat</th>\n",
       "      <th>Type</th>\n",
       "      <th>source</th>\n",
       "      <th>Presence</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>AMR</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Genomic Island</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Metal</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Novel Plasmid</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ERR1007500</td>\n",
       "      <td>United Kingdom</td>\n",
       "      <td>Agricultural</td>\n",
       "      <td>Type A</td>\n",
       "      <td>Phage</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Isolate       Geography       Habitat    Type          source  Presence\n",
       "0  ERR1007500  United Kingdom  Agricultural  Type A             AMR      18.0\n",
       "1  ERR1007500  United Kingdom  Agricultural  Type A  Genomic Island       4.0\n",
       "2  ERR1007500  United Kingdom  Agricultural  Type A           Metal       5.0\n",
       "3  ERR1007500  United Kingdom  Agricultural  Type A   Novel Plasmid       0.0\n",
       "4  ERR1007500  United Kingdom  Agricultural  Type A           Phage       3.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "intersect = counts[counts['Habitat'].isin(['Agricultural', 'Wastewater Mun.', 'Clinical'])]\n",
    "intersect.rename({'Country/Province': 'Geography'}, axis=1, inplace=True)\n",
    "intersect.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d700bf7c",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'model' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_5568/3470139634.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msummary\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'model' is not defined"
     ]
    }
   ],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2da58da0",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'source' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_5568/3189594945.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mintersect\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mintersect\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'source'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0msource\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mfactor_groups\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgroupby\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Type'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Habitat'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Geography'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'source' is not defined"
     ]
    }
   ],
   "source": [
    "t = intersect[intersect['source']==source]\n",
    "factor_groups = t.groupby(['Type', 'Habitat', 'Geography'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4a58f83e",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'source' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_5568/366932618.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msmf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mols\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Presence ~ C(Type) + C(Habitat) + C(Geography) + C(Type):C(Habitat) +C(Type):C(Geography) + C(Habitat):C(Geography) + C(Type):C(Geography):C(Habitat)'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mintersect\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mintersect\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'source'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0msource\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0md\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0manova_lm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtyp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0md\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'source' is not defined"
     ]
    }
   ],
   "source": [
    "model = smf.ols('Presence ~ C(Type) + C(Habitat) + C(Geography) + C(Type):C(Habitat) +C(Type):C(Geography) + C(Habitat):C(Geography) + C(Type):C(Geography):C(Habitat)', data=intersect[intersect['source']==source]).fit()\n",
    "d = sm.stats.anova_lm(model, typ=3)\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "c6734497",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "too many values to unpack (expected 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_20868/1769478609.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mfactor_groups\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Type'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Presence'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: too many values to unpack (expected 2)"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x432 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(6,6))\n",
    "for values, group in factor_groups:\n",
    "    i, j = values\n",
    "    plt.scatter(group['Type'], group['Presence'])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f81bd351",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AMR_p</th>\n",
       "      <th>Metal_p</th>\n",
       "      <th>VF_p</th>\n",
       "      <th>Plasmid_p</th>\n",
       "      <th>Genomic Island_p</th>\n",
       "      <th>Phage_p</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(Type)</th>\n",
       "      <td>6.423692e-04</td>\n",
       "      <td>0.236498</td>\n",
       "      <td>8.917313e-02</td>\n",
       "      <td>5.872894e-01</td>\n",
       "      <td>9.585507e-01</td>\n",
       "      <td>4.687334e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Habitat)</th>\n",
       "      <td>5.088880e-15</td>\n",
       "      <td>0.072621</td>\n",
       "      <td>1.288258e-09</td>\n",
       "      <td>1.736728e-31</td>\n",
       "      <td>5.027381e-40</td>\n",
       "      <td>1.142217e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Geography)</th>\n",
       "      <td>8.342733e-01</td>\n",
       "      <td>0.000042</td>\n",
       "      <td>2.277079e-04</td>\n",
       "      <td>8.756278e-03</td>\n",
       "      <td>6.459591e-01</td>\n",
       "      <td>7.798546e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Type):C(Habitat)</th>\n",
       "      <td>6.305390e-04</td>\n",
       "      <td>0.025376</td>\n",
       "      <td>1.462488e-04</td>\n",
       "      <td>1.054120e-04</td>\n",
       "      <td>5.263208e-07</td>\n",
       "      <td>7.414845e-04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Type):C(Geography)</th>\n",
       "      <td>3.035935e-01</td>\n",
       "      <td>0.689953</td>\n",
       "      <td>7.124783e-01</td>\n",
       "      <td>3.621372e-01</td>\n",
       "      <td>4.420687e-01</td>\n",
       "      <td>6.105836e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Habitat):C(Geography)</th>\n",
       "      <td>2.199025e-08</td>\n",
       "      <td>0.000162</td>\n",
       "      <td>6.576382e-01</td>\n",
       "      <td>2.504536e-04</td>\n",
       "      <td>1.763156e-13</td>\n",
       "      <td>1.163048e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Type):C(Geography):C(Habitat)</th>\n",
       "      <td>1.763488e-04</td>\n",
       "      <td>0.462426</td>\n",
       "      <td>4.881143e-03</td>\n",
       "      <td>1.743480e-01</td>\n",
       "      <td>2.612433e-03</td>\n",
       "      <td>7.036217e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                        AMR_p   Metal_p          VF_p  \\\n",
       "C(Type)                          6.423692e-04  0.236498  8.917313e-02   \n",
       "C(Habitat)                       5.088880e-15  0.072621  1.288258e-09   \n",
       "C(Geography)                     8.342733e-01  0.000042  2.277079e-04   \n",
       "C(Type):C(Habitat)               6.305390e-04  0.025376  1.462488e-04   \n",
       "C(Type):C(Geography)             3.035935e-01  0.689953  7.124783e-01   \n",
       "C(Habitat):C(Geography)          2.199025e-08  0.000162  6.576382e-01   \n",
       "C(Type):C(Geography):C(Habitat)  1.763488e-04  0.462426  4.881143e-03   \n",
       "\n",
       "                                    Plasmid_p  Genomic Island_p       Phage_p  \n",
       "C(Type)                          5.872894e-01      9.585507e-01  4.687334e-02  \n",
       "C(Habitat)                       1.736728e-31      5.027381e-40  1.142217e-08  \n",
       "C(Geography)                     8.756278e-03      6.459591e-01  7.798546e-01  \n",
       "C(Type):C(Habitat)               1.054120e-04      5.263208e-07  7.414845e-04  \n",
       "C(Type):C(Geography)             3.621372e-01      4.420687e-01  6.105836e-01  \n",
       "C(Habitat):C(Geography)          2.504536e-04      1.763156e-13  1.163048e-03  \n",
       "C(Type):C(Geography):C(Habitat)  1.743480e-01      2.612433e-03  7.036217e-02  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = []\n",
    "for source in intersect['source'].unique():\n",
    "    model = smf.ols('Presence ~ C(Type) + C(Habitat) + C(Geography) + C(Type):C(Habitat) +C(Type):C(Geography) + C(Habitat):C(Geography) + C(Type):C(Geography):C(Habitat)', data=intersect[intersect['source']==source]).fit()\n",
    "    d = sm.stats.anova_lm(model, typ=3)\n",
    "    d = d[['PR(>F)']]\n",
    "    d.rename({'PR(>F)': \"{}_p\".format(source)}, axis=1, inplace=True)\n",
    "    results.append(d)\n",
    "all_r = results[0]\n",
    "for d in results[1:]:\n",
    "    all_r = all_r.join(d)\n",
    "all_r.dropna().iloc[1:].dropna()[['AMR_p', 'Metal_p', 'VF_p', 'Plasmid_p', 'Genomic Island_p', 'Phage_p']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "14cbe0f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_r.dropna().iloc[1:].dropna()[['AMR_p', 'Metal_p', 'VF_p', 'Plasmid_p', 'Genomic Island_p', 'Phage_p']].to_csv('tables/Table1__3Way_ANOVA.csv', sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "14a5513a",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_r = all_r.dropna().iloc[1:].dropna()[['AMR_p', 'Metal_p', 'VF_p', 'Plasmid_p', 'Genomic Island_p', 'Phage_p']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78ac36c7",
   "metadata": {},
   "source": [
    "## Table 2: 2-Way ANOVA alberta data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1c04a7ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AMR_p</th>\n",
       "      <th>Metal_p</th>\n",
       "      <th>VF_p</th>\n",
       "      <th>Plasmid_p</th>\n",
       "      <th>Genomic Island_p</th>\n",
       "      <th>Phage_p</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(Habitat)</th>\n",
       "      <td>2.316110e-40</td>\n",
       "      <td>0.002480</td>\n",
       "      <td>3.392543e-14</td>\n",
       "      <td>2.194568e-62</td>\n",
       "      <td>1.610631e-55</td>\n",
       "      <td>3.953562e-11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Type)</th>\n",
       "      <td>2.337844e-07</td>\n",
       "      <td>0.120972</td>\n",
       "      <td>7.092580e-02</td>\n",
       "      <td>4.401716e-01</td>\n",
       "      <td>9.500725e-01</td>\n",
       "      <td>4.479566e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Habitat):C(Type)</th>\n",
       "      <td>4.711946e-08</td>\n",
       "      <td>0.000257</td>\n",
       "      <td>2.318734e-05</td>\n",
       "      <td>1.763849e-08</td>\n",
       "      <td>1.162846e-09</td>\n",
       "      <td>1.931045e-03</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           AMR_p   Metal_p          VF_p     Plasmid_p  \\\n",
       "C(Habitat)          2.316110e-40  0.002480  3.392543e-14  2.194568e-62   \n",
       "C(Type)             2.337844e-07  0.120972  7.092580e-02  4.401716e-01   \n",
       "C(Habitat):C(Type)  4.711946e-08  0.000257  2.318734e-05  1.763849e-08   \n",
       "\n",
       "                    Genomic Island_p       Phage_p  \n",
       "C(Habitat)              1.610631e-55  3.953562e-11  \n",
       "C(Type)                 9.500725e-01  4.479566e-02  \n",
       "C(Habitat):C(Type)      1.162846e-09  1.931045e-03  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts.rename({'Country/Province': \"Geography\"}, axis=1, inplace=True)\n",
    "results = []\n",
    "\n",
    "for source in counts['source'].unique():\n",
    "    model2 = smf.ols('Presence ~ C(Habitat) + C(Type) + C(Habitat):C(Type)', data=counts[(counts['source']==source) & (counts['Geography']=='Canada/Alberta')]).fit()\n",
    "    d = sm.stats.anova_lm(model2, typ=3)\n",
    "    d = d[['PR(>F)']]\n",
    "    d.rename({'PR(>F)': \"{}_p\".format(source)}, axis=1, inplace=True)\n",
    "    results.append(d)\n",
    "all_r2way = results[0]\n",
    "for d in results[1:]:\n",
    "    all_r2way = all_r2way.join(d)\n",
    "all_r2way.dropna().iloc[1:].dropna()[['AMR_p', 'Metal_p', 'VF_p', 'Plasmid_p', 'Genomic Island_p', 'Phage_p']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "08160a67",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_r2way = all_r2way.dropna().iloc[1:].dropna()[['AMR_p', 'Metal_p', 'VF_p', 'Plasmid_p', 'Genomic Island_p', 'Phage_p']]\n",
    "all_r2way.to_csv('tables/Table2__2Way_ANOVA.csv', sep=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "596b172d",
   "metadata": {},
   "source": [
    "## Table X (unassigned): Bonferroni Corrected t-test post-hocs\n",
    "My assumed procedure:\n",
    "1. Do the ANOVA\n",
    "2. For cells in the ANOVA that have significant differences, we can do a t-test post-hoc.\n",
    "\n",
    "Since the 2-way anova didn't have any insig. cells that were sig in three-way, just use the three-way.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "22dbcf1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "07ba6de8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POST-HOC FOR AMR, Type\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_5568/173186406.py:22: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['groups'] = groups\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POST-HOC FOR AMR, Habitat\n",
      "POST-HOC FOR AMR, Type_Habitat\n",
      "POST-HOC FOR AMR, Habitat_Geography\n",
      "POST-HOC FOR AMR, Type_Geography_Habitat\n",
      "POST-HOC FOR Metal, Geography\n",
      "POST-HOC FOR Metal, Type_Habitat\n",
      "POST-HOC FOR Metal, Habitat_Geography\n",
      "POST-HOC FOR VF, Habitat\n",
      "POST-HOC FOR VF, Geography\n",
      "POST-HOC FOR VF, Type_Habitat\n",
      "POST-HOC FOR VF, Type_Geography_Habitat\n",
      "POST-HOC FOR Plasmid, Habitat\n",
      "POST-HOC FOR Plasmid, Geography\n",
      "POST-HOC FOR Plasmid, Type_Habitat\n",
      "POST-HOC FOR Plasmid, Habitat_Geography\n",
      "POST-HOC FOR Genomic Island, Habitat\n",
      "POST-HOC FOR Genomic Island, Type_Habitat\n",
      "POST-HOC FOR Genomic Island, Habitat_Geography\n",
      "POST-HOC FOR Genomic Island, Type_Geography_Habitat\n",
      "POST-HOC FOR Phage, Type\n",
      "POST-HOC FOR Phage, Habitat\n",
      "POST-HOC FOR Phage, Type_Habitat\n",
      "POST-HOC FOR Phage, Habitat_Geography\n"
     ]
    }
   ],
   "source": [
    "factors = [['Type'],\n",
    "           ['Habitat'],\n",
    "           ['Geography'],\n",
    "           ['Type', 'Habitat'], \n",
    "           ['Type', 'Geography'],\n",
    "           ['Habitat', 'Geography'],\n",
    "           ['Type', 'Geography', 'Habitat']]\n",
    "\n",
    "def query(factor):\n",
    "    return ':'.join(['C({})'.format(f) for f in factor])\n",
    "\n",
    "tables3way = {}\n",
    "#First, for 3-way\n",
    "for source in ['AMR', 'Metal', 'VF', 'Plasmid', 'Genomic Island', 'Phage']:\n",
    "    data = intersect[intersect['source']==source]\n",
    "    for factor in factors:\n",
    "        if all_r.loc[query(factor), '{}_p'.format(source)] <= 0.05:\n",
    "            print(\"POST-HOC FOR {0}, {1}\".format(source, '_'.join(factor)))\n",
    "            groups = []\n",
    "            for i, row in data.iterrows():\n",
    "                groups.append(' X '.join(row[factor].values))\n",
    "            data['groups'] = groups\n",
    "            tu = sp.posthoc_tukey(data[['groups', 'Presence']], group_col='groups', val_col='Presence')\n",
    "            \"\"\"\n",
    "            df = pd.DataFrame(tu.summary())\n",
    "            header = df.iloc[0]\n",
    "            header = [str(h) for h in header]\n",
    "            df = df[1:]\n",
    "            df.columns = header\n",
    "            \"\"\"\n",
    "            tables3way['{0}_{1}'.format(source, '_'.join(factor))] = tu\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "21ec87cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, table in tables.items():\n",
    "    table.to_csv('tables/posthoc_tests/{}_3way_post.csv'.format(key), sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "000dccdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AMR_p</th>\n",
       "      <th>Metal_p</th>\n",
       "      <th>VF_p</th>\n",
       "      <th>Plasmid_p</th>\n",
       "      <th>Genomic Island_p</th>\n",
       "      <th>Phage_p</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(Habitat)</th>\n",
       "      <td>2.316110e-40</td>\n",
       "      <td>0.002480</td>\n",
       "      <td>3.392543e-14</td>\n",
       "      <td>2.194568e-62</td>\n",
       "      <td>1.610631e-55</td>\n",
       "      <td>3.953562e-11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Type)</th>\n",
       "      <td>2.337844e-07</td>\n",
       "      <td>0.120972</td>\n",
       "      <td>7.092580e-02</td>\n",
       "      <td>4.401716e-01</td>\n",
       "      <td>9.500725e-01</td>\n",
       "      <td>4.479566e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(Habitat):C(Type)</th>\n",
       "      <td>4.711946e-08</td>\n",
       "      <td>0.000257</td>\n",
       "      <td>2.318734e-05</td>\n",
       "      <td>1.763849e-08</td>\n",
       "      <td>1.162846e-09</td>\n",
       "      <td>1.931045e-03</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           AMR_p   Metal_p          VF_p     Plasmid_p  \\\n",
       "C(Habitat)          2.316110e-40  0.002480  3.392543e-14  2.194568e-62   \n",
       "C(Type)             2.337844e-07  0.120972  7.092580e-02  4.401716e-01   \n",
       "C(Habitat):C(Type)  4.711946e-08  0.000257  2.318734e-05  1.763849e-08   \n",
       "\n",
       "                    Genomic Island_p       Phage_p  \n",
       "C(Habitat)              1.610631e-55  3.953562e-11  \n",
       "C(Type)                 9.500725e-01  4.479566e-02  \n",
       "C(Habitat):C(Type)      1.162846e-09  1.931045e-03  "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_r2way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "83c71fa9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POST-HOC FOR AMR, Type\n",
      "POST-HOC FOR AMR, Habitat\n",
      "POST-HOC FOR AMR, Habitat_Type\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_5568/387158977.py:15: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['groups'] = groups\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POST-HOC FOR Metal, Habitat\n",
      "POST-HOC FOR Metal, Habitat_Type\n",
      "POST-HOC FOR VF, Habitat\n",
      "POST-HOC FOR VF, Habitat_Type\n",
      "POST-HOC FOR Plasmid, Habitat\n",
      "POST-HOC FOR Plasmid, Habitat_Type\n",
      "POST-HOC FOR Genomic Island, Habitat\n",
      "POST-HOC FOR Genomic Island, Habitat_Type\n",
      "POST-HOC FOR Phage, Type\n",
      "POST-HOC FOR Phage, Habitat\n",
      "POST-HOC FOR Phage, Habitat_Type\n"
     ]
    }
   ],
   "source": [
    "factors = [['Type'],\n",
    "           ['Habitat'],\n",
    "           ['Habitat','Type'],]\n",
    "\n",
    "tables2way = {}\n",
    "#repeat for 2-way\n",
    "for source in ['AMR', 'Metal', 'VF', 'Plasmid', 'Genomic Island', 'Phage']:\n",
    "    data=counts[(counts['source']==source) & (counts['Geography']=='Canada/Alberta')]\n",
    "    for factor in factors:\n",
    "        if all_r2way.loc[query(factor), '{}_p'.format(source)] <= 0.05:\n",
    "            print(\"POST-HOC FOR {0}, {1}\".format(source, '_'.join(factor)))\n",
    "            groups = []\n",
    "            for i, row in data.iterrows():\n",
    "                groups.append(' X '.join(row[factor].values))\n",
    "            data['groups'] = groups\n",
    "            tu = sp.posthoc_tukey(data[['groups', 'Presence']], group_col='groups', val_col='Presence')\n",
    "            \"\"\"\n",
    "            df = pd.DataFrame(tu.summary())\n",
    "            header = df.iloc[0]\n",
    "            header = [str(h) for h in header]\n",
    "            df = df[1:]\n",
    "            df.columns = header\n",
    "            \"\"\"\n",
    "            tables2way['{0}_{1}'.format(source, '_'.join(factor))] = tu\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "7fe3e10d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, table in tables.items():\n",
    "    table.to_csv('tables/posthoc_tests/{}_2way_post.csv'.format(key), sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "25f00cf4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['AMR_Type', 'AMR_Habitat', 'AMR_Habitat_Type', 'Metal_Habitat', 'Metal_Habitat_Type', 'VF_Habitat', 'VF_Habitat_Type', 'Plasmid_Habitat', 'Plasmid_Habitat_Type', 'Genomic Island_Habitat', 'Genomic Island_Habitat_Type', 'Phage_Type', 'Phage_Habitat', 'Phage_Habitat_Type'])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tables2way.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "ba829b80",
   "metadata": {},
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('tables/posthoc_spreadsheet.xlsx',engine='xlsxwriter')\n",
    "workbook = writer.book\n",
    "\n",
    "significant = workbook.add_format({'bg_color': '#93c47d'})\n",
    "medium = workbook.add_format({'bg_color': '#ffd966'})\n",
    "bad = workbook.add_format({'bg_color': \"#e06666\"})\n",
    "\n",
    "for source in ['AMR', 'Metal', 'VF', 'Plasmid', 'Genomic Island', 'Phage']:\n",
    "    offset = 0\n",
    "    worksheet = workbook.add_worksheet(source)\n",
    "    writer.sheets[source] = worksheet\n",
    "    for key in tables3way.keys():\n",
    "        if key.startswith(source):\n",
    "            worksheet.write_string(offset,0, '{0} {1} {2}'.format(source, key.replace('_',' X '), \"3-way ANOVOA posthoc, Tukey's HSD\"))\n",
    "            offset += 1\n",
    "            tables3way[key].to_excel(writer, sheet_name=source, startrow=offset, startcol=0)\n",
    "            offset += tables3way[key].shape[0] + 4\n",
    "    for key in tables2way.keys():\n",
    "        if key.startswith(source):\n",
    "            worksheet.write_string(offset,0, '{0} {1} {2}'.format(source, key.replace('_', ' X '), \"2-way ANOVOA posthoc, Tukey's HSD\"))\n",
    "            offset += 1\n",
    "            tables2way[key].to_excel(writer, sheet_name=source, startrow=offset, startcol=0)\n",
    "            offset += tables2way[key].shape[0] + 4\n",
    "    worksheet.conditional_format('A1:Z1000', {'type': 'cell', \n",
    "                                              'criteria': '<=',\n",
    "                                              'value': 0.05,\n",
    "                                              'format': significant})\n",
    "    worksheet.conditional_format('A1:Z1000', {'type': 'cell', \n",
    "                                          'criteria': 'between',\n",
    "                                          'minimum': 0.05,\n",
    "                                          'maximum': 0.15,\n",
    "                                          'format': medium})\n",
    "    worksheet.conditional_format('A1:Z1000', {'type': 'cell',\n",
    "                                              'criteria': '>',\n",
    "                                              'value': 0.15,\n",
    "                                              'format':bad})\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc32266d",
   "metadata": {},
   "source": [
    "## Table S3 Presence/absence of everything"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "3037ff84",
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_gene(percent):\n",
    "    if percent >= 0.99:\n",
    "        return \"Core\"\n",
    "    elif percent >= 0.95:\n",
    "        return \"Soft Core\"\n",
    "    elif percent >= 0.15:\n",
    "        return \"Shell\"\n",
    "    else:\n",
    "        return \"Cloud\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "45ac9001",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['AMR', 'Genomic Island', 'Metal', 'Novel Plasmid', 'Phage',\n",
       "       'Plasmid', 'VF'], dtype=object)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts['source'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "6d4f89a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "frames =[]\n",
    "for feat in ['AMR', 'Metal', 'VF', 'Plasmid', 'Genomic Island', 'Phage']:\n",
    "    table = pd.read_table('tables/feature_presence_percentages/AMR_summary.csv', sep=',')[['Feature', 'Presence', 'Percentage']]\n",
    "    table['source'] = feat\n",
    "    frames.append(table)\n",
    "type_counts = pd.concat(frames)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "75e0743f",
   "metadata": {},
   "outputs": [],
   "source": [
    "type_counts['gene_type'] = type_counts['Percentage'].map(assign_gene)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "0c3f1be4",
   "metadata": {},
   "outputs": [],
   "source": [
    "type_counts.set_index('Feature').to_csv('tables/TableX__target_feature_core_assignment.csv', sep=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "6bacbc07",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the roary gene presence absence table to count genes.\n",
    "gpa = pd.read_table('gene_presence_absence_roary.csv', sep = ',', low_memory=False)\n",
    "genomes = gpa.columns[14:] #14th column is start of genomes\n",
    "# Convert to presence absence, remove reference and outgroup\n",
    "roary = gpa.set_index('Gene')[genomes].T.notnull().astype(int).loc[meta['Isolate']]\n",
    "\n",
    "roary.columns.name = ''\n",
    "roary.index.name = 'Isolate'\n",
    "\n",
    "#convert P/A to long form\n",
    "roary_long = pd.melt(roary.reset_index(), id_vars=['Isolate'], var_name='Feature', value_name='Presence')\n",
    "\n",
    "# count genomes\n",
    "roary_t = roary_long.groupby('Feature').sum()\n",
    "# percentage of genomes\n",
    "roary_t = roary_t.join(roary_t.div(1273).rename({'Presence': 'Percentage'}, axis=1))\n",
    "roary_t.sort_values(by='Presence', ascending=False, inplace=True)\n",
    "# Annotate with core/cloud/shell\n",
    "roary_t['source'] = 'Roary'\n",
    "roary_t['gene_type'] = roary_t['Percentage'].map(assign_gene)\n",
    "\n",
    "roary_t.to_csv('tables/TableX__roary_feature_core_assignment.csv', sep=',')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "18e0ee75",
   "metadata": {},
   "outputs": [],
   "source": [
    "pa = pd.read_table('all_features_presence_absence.tsv', sep='\\t', index_col=0)\n",
    "all_pa = pa.join(roary, rsuffix='_roary')\n",
    "all_pa.to_csv('tables/TableS3__all_presence_absence.csv', sep=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ad40964",
   "metadata": {},
   "source": [
    "## Table S5 (?): Pagel vs labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "id": "39399443",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Isolate</th>\n",
       "      <th>Feature</th>\n",
       "      <th>Presence</th>\n",
       "      <th>Habitat</th>\n",
       "      <th>Country/Province</th>\n",
       "      <th>Clade</th>\n",
       "      <th>source</th>\n",
       "      <th>Sampling Location</th>\n",
       "      <th>Geography</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SRR14026455</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Wastewater Agr.</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade B</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Wastewater Agr. Ab</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SRR14026462</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SRR14026496</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SRR14026467</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade A</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SRR14026460</td>\n",
       "      <td>AA083</td>\n",
       "      <td>0.0</td>\n",
       "      <td>Natural Water</td>\n",
       "      <td>Canada/Alberta</td>\n",
       "      <td>Clade B</td>\n",
       "      <td>Plasmid</td>\n",
       "      <td>Natural Water Ab</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Isolate Feature  Presence          Habitat Country/Province    Clade  \\\n",
       "0  SRR14026455   AA083       0.0  Wastewater Agr.   Canada/Alberta  Clade B   \n",
       "1  SRR14026462   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "2  SRR14026496   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "3  SRR14026467   AA083       0.0    Natural Water   Canada/Alberta  Clade A   \n",
       "4  SRR14026460   AA083       0.0    Natural Water   Canada/Alberta  Clade B   \n",
       "\n",
       "    source   Sampling Location Geography  \n",
       "0  Plasmid  Wastewater Agr. Ab        Ab  \n",
       "1  Plasmid    Natural Water Ab        Ab  \n",
       "2  Plasmid    Natural Water Ab        Ab  \n",
       "3  Plasmid    Natural Water Ab        Ab  \n",
       "4  Plasmid    Natural Water Ab        Ab  "
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fe = pd.read_table('genomes_vs_features_longform.csv', sep=',')\n",
    "loc_abbr_map = {'United Kingdom': 'UK', \n",
    "                'Canada/Alberta': 'Ab'}\n",
    "\n",
    "fe['Sampling Location'] = fe.apply(lambda row: ' '.join([row['Habitat'], loc_abbr_map[row['Country/Province']]]), axis=1)\n",
    "fe['Geography']= [loc_abbr_map[c] for c in fe['Country/Province']]\n",
    "fe.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "id": "fd8905a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rename_clade(string):\n",
    "    if string == 'CladeB':\n",
    "        return 'Type B'\n",
    "    if type(string)==str:\n",
    "        return string.replace('Clade', 'Type')\n",
    "    return string\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "id": "733d6258",
   "metadata": {},
   "outputs": [],
   "source": [
    "fe['Type'] = fe['Clade'].map(rename_clade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "id": "2528b0ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Metadata Type</th>\n",
       "      <th>Metadata</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Habitat</td>\n",
       "      <td>Wastewater Agr.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Habitat</td>\n",
       "      <td>Natural Water</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Habitat</td>\n",
       "      <td>Wastewater Mun.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Habitat</td>\n",
       "      <td>Clinical</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Habitat</td>\n",
       "      <td>Agricultural</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Geography</td>\n",
       "      <td>Ab</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Geography</td>\n",
       "      <td>UK</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Type</td>\n",
       "      <td>Type B</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Type</td>\n",
       "      <td>Type A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Type</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Metadata Type         Metadata\n",
       "0       Habitat  Wastewater Agr.\n",
       "1       Habitat    Natural Water\n",
       "2       Habitat  Wastewater Mun.\n",
       "3       Habitat         Clinical\n",
       "4       Habitat     Agricultural\n",
       "5     Geography               Ab\n",
       "6     Geography               UK\n",
       "7          Type           Type B\n",
       "8          Type           Type A\n",
       "9          Type              NaN"
      ]
     },
     "execution_count": 189,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To sort properly \n",
    "recs = []\n",
    "for col in ['Habitat', 'Geography', 'Type']:\n",
    "    for value in fe[col].unique():\n",
    "        recs.append({\"Metadata Type\": col,\n",
    "                     \"Metadata\": value})\n",
    "        \n",
    "legend = pd.DataFrame.from_records(recs)\n",
    "legend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "id": "43c9f2e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "l1 = ['Agriculture',\n",
    " 'Clinical',\n",
    " 'Natural Water',\n",
    " 'Wastewater Agr.',\n",
    " 'Wastewater Municipal']\n",
    "l2 = ['Agricultural',\n",
    " 'Clinical',\n",
    " 'Natural Water',\n",
    " 'Wastewater Agr.',\n",
    " 'Wastewater Mun.']\n",
    "\n",
    "rename = {k:v for k, v in zip(l1,l2)}\n",
    "rename['Clade A'] = 'Type A'\n",
    "rename['Clade B'] = 'Type B'\n",
    "rename['CladeB'] = 'Type B'\n",
    "rename = {**rename, ** loc_abbr_map, \"Alberta\": \"Ab\"}\n",
    "def fix_pagel(val):\n",
    "    if val in rename.keys():\n",
    "        return rename[val]\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "8bf358d7",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'fe' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_20868/3428163102.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mlr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_table\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'pagel_vs_metadata_LR.csv'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m','\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mfeats\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfe\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Feature'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'source'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdrop_duplicates\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'fe' is not defined"
     ]
    }
   ],
   "source": [
    "p_val = pd.read_table('pagel_vs_metadata_pvalue.csv', sep=',')\n",
    "lr = pd.read_table('pagel_vs_metadata_LR.csv', sep=',')\n",
    "\n",
    "feats = fe[['Feature', 'source']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "12512e76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.11022302462516e-16"
      ]
     },
     "execution_count": 172,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_table('../../Indizio/data/pagel_features_p.csv', sep=',', index_col=0)\n",
    "v = df.values\n",
    "np.nanmin(v[np.nonzero(v)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "id": "aa8bff80",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_val['Unnamed: 0'] = p_val['Unnamed: 0'].map(fix_pagel)\n",
    "lr['Unnamed: 0'] = lr['Unnamed: 0'].map(fix_pagel)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "id": "dfcaf347",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "for feat in feats['source'].unique():\n",
    "    f = feats[feats['source']==feat]\n",
    "    f_p = p_val[['Unnamed: 0'] + list(f[f['Feature'].isin(p_val.columns)]['Feature'])]\n",
    "    f_lr = lr[['Unnamed: 0'] + list(f[f['Feature'].isin(lr.columns)]['Feature'])]\n",
    "    \n",
    "    f_pg = f_p.melt(id_vars=['Unnamed: 0'])\n",
    "    f_pg.rename({'Unnamed: 0': 'Metadata', 'variable': 'Feature', 'value': 'p value'}, axis=1, inplace=True)\n",
    "    f_pg.set_index(['Feature', 'Metadata'], inplace=True)\n",
    "    \n",
    "    f_lr = lr[['Unnamed: 0'] + list(f[f['Feature'].isin(lr.columns)]['Feature'])]\n",
    "    f_lrg = f_lr.melt(id_vars=['Unnamed: 0'])\n",
    "    f_lrg.rename({'Unnamed: 0': 'Metadata', 'variable': 'Feature', 'value': 'LR'}, axis=1, inplace=True)\n",
    "    f_lrg.set_index(['Feature', 'Metadata'], inplace=True)\n",
    "    f_pg.join(f_lrg).to_csv('tables/pagel_vs_habitat/{}_vs_habitat_pagel.csv'.format(feat), sep=',')\n",
    "    f_pg_lr = f_pg.join(f_lrg)\n",
    "    data = fe[fe['source']==feat]\n",
    "    perc_frames = []\n",
    "    presence_frames = []\n",
    "\n",
    "    total = data.groupby('Feature').sum()\n",
    "    total = total.join(total.div(1273).rename({'Presence': \"Percentage\"}, axis=1))\n",
    "    total.sort_values(by='Presence', ascending=False, inplace=True)\n",
    "\n",
    "    for col in ['Type', 'Habitat', 'Geography']:\n",
    "        p = data.groupby([col, 'Feature']).agg({'Presence': 'sum'})\n",
    "        c = data.groupby([col]).agg({'Isolate': 'nunique'})\n",
    "        c.rename({'Isolate': 'Presence'}, axis=1, inplace=True)\n",
    "        perc = p.div(c, level=col).reset_index().pivot(index='Feature', columns=col, values='Presence')\n",
    "        presence = p.reset_index().pivot(index='Feature', columns=col, values='Presence')\n",
    "\n",
    "        perc.columns.name = ''\n",
    "        presence.columns.name =''\n",
    "\n",
    "        perc_frames.append(perc.reset_index().melt(id_vars=['Feature']).rename({'': 'Metadata', 'value': 'Proportion Present'},axis=1).set_index(['Feature', 'Metadata']))\n",
    "        presence_frames.append(presence.reset_index().melt(id_vars=['Feature']).rename({'': 'Metadata', 'value': 'n Genomes Present'},axis=1).set_index(['Feature', 'Metadata']))\n",
    "\n",
    "    final = f_pg_lr.reset_index().join(legend.set_index('Metadata'), on ='Metadata')\\\n",
    "        .sort_values(by=['Feature', 'Metadata Type', 'Metadata'])\\\n",
    "        .drop('Metadata Type',axis=1)\\\n",
    "        .set_index(['Feature', 'Metadata'])\\\n",
    "        .join(pd.concat(perc_frames)).join(pd.concat(presence_frames))\n",
    "    final.to_csv('tables/pagel_vs_habitat/{}_pagel_with_counts.csv'.format(feat), sep=',')\n",
    "    results[feat] = final"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}