rebuilt_variance_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "%load_ext rpy2.ipython\n",
    "\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import warnings\n",
    "import numpy as np\n",
    "import microbe_directory as md\n",
    "import seaborn as sns\n",
    "\n",
    "from matplotlib import pyplot as plt        \n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "\n",
    "from capalyzer.packet_parser import DataTableFactory, NCBITaxaTree, annotate_taxa, TaxaTree\n",
    "from capalyzer.packet_parser.data_utils import group_small_cols\n",
    "from capalyzer.packet_parser.diversity_metrics import (\n",
    "    shannon_entropy, richness, chao1, rarefaction_analysis\n",
    ")\n",
    "from sklearn.decomposition import PCA\n",
    "from scipy.cluster.hierarchy import linkage, cophenet, leaves_list\n",
    "from scipy.spatial.distance import squareform, pdist, jensenshannon\n",
    "\n",
    "from os.path import join\n",
    "from metasub_utils.packet_parse import MetaSUBTableFactory\n",
    "from metasub_utils.packet_parse.metadata_ontology import add_ontology\n",
    "from capalyzer.packet_parser.experimental import umap\n",
    "from capalyzer.packet_parser.data_utils import group_small_cols\n",
    "from capalyzer.packet_parser.normalize import proportions, prevalence\n",
    "from plotnine import *\n",
    "from scipy.cluster.hierarchy import fcluster\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from capalyzer.constants import MICROBE_DIR\n",
    "\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import svm \n",
    "from sklearn.metrics import (\n",
    "    precision_score,\n",
    "    recall_score,\n",
    "    confusion_matrix,\n",
    "    classification_report,\n",
    "    accuracy_score,\n",
    ")\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "NAN = float('nan')\n",
    "PACKET_DIR = '/home/dcdanko/Dropbox/resources_and_shared/metasub_data_packets/release_packet_jun12_2019/'\n",
    "PACKET_DIR = '/home/dcdanko/Dropbox/resources_and_shared/metasub_data_packets/data_packet_2020_02_23_16_12'\n",
    "ncbi_tree = NCBITaxaTree.parse_files()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_title(el):\n",
    "    sel = str(el)\n",
    "    if not el:\n",
    "        return el\n",
    "    tkns = sel.split('_')\n",
    "    tkns = [tkn[0].upper() + tkn[1:] for tkn in tkns if tkn]\n",
    "    return ' '.join(tkns)\n",
    "\n",
    "\n",
    "def categorical_continents(continents):\n",
    "    return pd.Categorical(\n",
    "        pd.Series(continents).map(to_title),\n",
    "        categories=[\n",
    "            'North America',\n",
    "            'East Asia',\n",
    "            'Europe',\n",
    "            'Sub Saharan Africa',\n",
    "            'South America',\n",
    "            'Middle East',\n",
    "            'Oceania',\n",
    "            'Nan',\n",
    "        ],\n",
    "        ordered=True,\n",
    "    )\n",
    "\n",
    "def normalize_surface(el):\n",
    "    try:\n",
    "        el = el.lower()\n",
    "    except AttributeError:\n",
    "        return el\n",
    "    el = '_'.join(el.split())\n",
    "    return {\n",
    "         'ground': 'floor',\n",
    "         'palm_left': 'human_hand',\n",
    "         'vert_pole': 'pole',\n",
    "         'palm_right': 'human_hand',\n",
    "         'pedestrian_crossing_button': 'button',\n",
    "         'wooden bench': 'seat',\n",
    "         'gargabe_can': 'garbage',\n",
    "         'horiz_pole': 'railing',\n",
    "         'stairwell railing': 'railing',\n",
    "         'bike kiosk': 'kiosk',\n",
    "         'bench1': 'seat',\n",
    "         'kiosk2': 'kiosk',\n",
    "         'kiosk1': 'kiosk',\n",
    "         'ticket_machine': 'kiosk',\n",
    "         'escalator_handrail': 'railing',\n",
    "         'bench2': 'seat',\n",
    "         'lift_buttons': 'button',\n",
    "         'overhead_handrail': 'railing',\n",
    "         'poll': 'pole',\n",
    "        'garbage_can': 'garbage',\n",
    "        'bench': 'seat',\n",
    "        'door_handle;door_handle': 'door_knob',\n",
    "        'glass': 'window',\n",
    "        'ceiling_rail': 'railing',\n",
    "        'seat_rail': 'railing',\n",
    "        'ticket_hall;ticket_machine': 'kiosk',\n",
    "        'bench;platform': 'seat',\n",
    "        'platform;bench': 'seat',\n",
    "        'ticket_machine;ticket_hall': 'kiosk',\n",
    "        'wooden_bench': 'seat',\n",
    "        'stairwell_railing': 'railing',\n",
    "        'trolley_handle': '',\n",
    "        'turnstile_or_alternatives;turnstile_or_alternatives': 'turnstile',\n",
    "        'ticket_kiosks;ticket_kiosks': 'kiosk',\n",
    "        'ticket_machine;ticket_machine': 'kiosk',\n",
    "        'station_eletronic_kiosk': 'kiosk',\n",
    "        'station_railing': 'railing',\n",
    "        'station_seat': 'seat',\n",
    "        'center_seat': 'seat',\n",
    "        'negative_control_(air)': '',\n",
    "        'ceiling_rail;ceiling_rail': 'railing',\n",
    "        'rail': 'railing',\n",
    "        'seat_near_door': 'seat',\n",
    "        'ticketing_machine;ticketing_machine_': 'kiosk',\n",
    "        'jetway_1_seats_waiting_area': 'seat',\n",
    "    }.get(el, el)\n",
    "\n",
    "def add_ontology(metadata):\n",
    "    \"\"\"Return a pandas dataframe with metadata and MetaSUB ontologies.\"\"\"\n",
    "    metadata = add_surface_ontology(metadata)\n",
    "    metadata = add_place_ontology(metadata)\n",
    "    return metadata\n",
    "\n",
    "\n",
    "def has_keyword(target, *queries):\n",
    "    try:\n",
    "        for query in queries:\n",
    "            if query in target:\n",
    "                return True\n",
    "    except TypeError:\n",
    "        pass\n",
    "    return False\n",
    "\n",
    "\n",
    "def add_surface_ontology(metadata):\n",
    "    \"\"\"Return a pandas dataframe with metadata and surface ontologies.\"\"\"\n",
    "    metadata, tbl = metadata.copy(), {}\n",
    "    for val in metadata['surface_material'].unique():\n",
    "        if has_keyword(val, 'glass', 'metal', 'steel', 'copper'):\n",
    "            tbl[val] = ('metal', 'impermeable')\n",
    "        elif has_keyword(val, 'stone', 'marble', 'ceramic', 'concrete', 'cement', 'granite'):\n",
    "            tbl[val] = ('stone', 'impermeable')\n",
    "        elif has_keyword(val, 'plastic', 'rubber', 'vinyl', 'pvc', 'formica'):\n",
    "            tbl[val] = ('plastic', 'impermeable')\n",
    "        elif has_keyword(val, 'fabric', 'cloth', 'carpet'):\n",
    "            tbl[val] = ('fabric', 'permeable')\n",
    "        elif has_keyword(val, 'hand', 'flesh', 'wood', 'leather', 'fiber'):\n",
    "            tbl[val] = ('biological', 'permeable')\n",
    "        elif has_keyword(val, 'control'):\n",
    "            tbl[val] = ('control', 'control')\n",
    "        else:\n",
    "            tbl[val] = (NAN, NAN)\n",
    "    metadata['surface_ontology_fine'] = metadata['surface_material'].apply(lambda x: tbl[x][0])\n",
    "    metadata['surface_ontology_coarse'] = metadata['surface_material'].apply(lambda x: tbl[x][1])\n",
    "    return metadata\n",
    "\n",
    "def add_place_ontology(metadata):\n",
    "    \"\"\"Return a pandas dataframe with metadata and place ontologies.\"\"\"\n",
    "    metadata = metadata.copy()\n",
    "    metadata['coastal'] = metadata.apply(coastal, result_type='reduce', axis=1)\n",
    "    return metadata\n",
    "\n",
    "\n",
    "def coastal(row):\n",
    "    if row['coastal_city'] == 'yes':\n",
    "        return ('coastal', 'coastal')\n",
    "    if float(row['city_elevation_meters']) > 1000:\n",
    "        return ('high altitude', 'not_coastal')\n",
    "    return ('low altitude', 'not_coastal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5132, 41)\n"
     ]
    }
   ],
   "source": [
    "class MetaSUBFiguresData:\n",
    "\n",
    "    def __init__(self, packet_dir, ncbi_tree=None):\n",
    "        self.tabler = DataTableFactory(packet_dir, metadata_tbl='metadata/complete_metadata.csv')\n",
    "        print(self.tabler.metadata.shape)\n",
    "        self.tabler.metadata = add_ontology(self.tabler.metadata)\n",
    "        self.tabler.metadata['city_elevation'] = self.tabler.metadata['coastal'].map(lambda el: el[1])\n",
    "        self.tabler.metadata['coastal'] = self.tabler.metadata['coastal'].map(lambda el: el[0])\n",
    "        self.tabler.metadata['continent'] = categorical_continents(self.tabler.metadata['continent'])\n",
    "        self.tabler.metadata['city'] = self.tabler.metadata['city'].map(to_title)\n",
    "        no_air = self.tabler.metadata.query('project != \"CSD17_AIR\"')\n",
    "        self.tabler = self.tabler.copy(new_metadata=no_air)\n",
    "        self.meta = self.tabler.metadata\n",
    "        self._long_taxa = None\n",
    "        self._wide_taxa = None\n",
    "\n",
    "    @property\n",
    "    def amr_genes(self):\n",
    "        return self.build_amr_genes()\n",
    "\n",
    "    @property\n",
    "    def long_taxa(self):\n",
    "        if self._long_taxa is None:\n",
    "            self._long_taxa = self.build_long_taxa()\n",
    "        return self._long_taxa\n",
    "\n",
    "    @property\n",
    "    def wide_taxa(self):\n",
    "        if self._wide_taxa is None:\n",
    "            self._wide_taxa = self.build_wide_taxonomy()\n",
    "        return self._wide_taxa\n",
    "\n",
    "    @property\n",
    "    def wide_taxa_rel(self):\n",
    "        return proportions(self.wide_taxa)\n",
    "    \n",
    "    def build_amr_genes(self):\n",
    "        \"\"\"TODO: REVIEW.\"\"\"\n",
    "        amrs = self.tabler.amrs(kind='gene', remove_zero_rows=False).dropna()\n",
    "        amrs = (amrs.T / (amrs.T.sum() + 0.000001)).T\n",
    "        return amrs\n",
    "        \n",
    "    def build_wide_taxonomy(self):\n",
    "        \"\"\"Return a pandas df with species in columns, samples in rows. Values are read counts.\"\"\"\n",
    "        # return self.tabler.taxonomy()\n",
    "        return self.long_taxa.pivot(index='sample', columns='taxa_name', values='reads').fillna(0)\n",
    "\n",
    "    def build_long_taxa(self):\n",
    "        long_taxa = pd.read_csv(join(self.tabler.packet_dir, 'taxonomy/refseq.krakenhll_longform.csv.gz'))\n",
    "        long_taxa = long_taxa.rename(columns={\n",
    "            'Unnamed: 0': 'sample',\n",
    "            'Unnamed: 1': 'taxa_name',\n",
    "            'Unnamed: 2': 'taxa_id',\n",
    "            'Unnamed: 3': 'rank',\n",
    "            'Unnamed: 4': 'mpa',\n",
    "        })\n",
    "        sample_list = self.tabler.metadata.index\n",
    "        long_taxa = long_taxa.loc[long_taxa['sample'].isin(sample_list)]\n",
    "        long_taxa = long_taxa.query('rank == \"species\"')\n",
    "        for col in ['reads', 'kmers', 'tax_reads']:\n",
    "            long_taxa[col] = long_taxa[col].fillna(0).map(int)\n",
    "        for col in ['dup', 'cov', 'percent']:\n",
    "            long_taxa[col] = long_taxa[col].fillna(0).map(float)\n",
    "\n",
    "        long_taxa['kmer_read_ratio'] = long_taxa['kmers'] / long_taxa['reads']\n",
    "\n",
    "        long_taxa = long_taxa.query('kmer_read_ratio >= 2.5')\n",
    "        long_taxa = long_taxa.query('kmers >= 64')\n",
    "        long_taxa = long_taxa.query('reads >= 3')\n",
    "        long_taxa = long_taxa.query('taxa_name != \"Human endogenous retrovirus\"')\n",
    "        long_taxa = long_taxa.query('taxa_name != \"Human endogenous retrovirus K\"')\n",
    "        long_taxa = long_taxa.loc[long_taxa['sample'].isin(self.meta.index)]\n",
    "        return long_taxa\n",
    "        \n",
    "data = MetaSUBFiguresData(PACKET_DIR)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>metasub_name</th>\n",
       "      <th>core_project</th>\n",
       "      <th>project</th>\n",
       "      <th>city</th>\n",
       "      <th>city_code</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>surface_material</th>\n",
       "      <th>control_type</th>\n",
       "      <th>elevation</th>\n",
       "      <th>...</th>\n",
       "      <th>index_sequence</th>\n",
       "      <th>location_type</th>\n",
       "      <th>hudson_alpha_uid</th>\n",
       "      <th>other_project_uid</th>\n",
       "      <th>sample_type</th>\n",
       "      <th>sl_name</th>\n",
       "      <th>surface_ontology_fine</th>\n",
       "      <th>surface_ontology_coarse</th>\n",
       "      <th>coastal</th>\n",
       "      <th>city_elevation</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>uuid</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>haib17KIU4866_H7HJMCCXY_SL272833</th>\n",
       "      <td>CSD16-LIS-CTRL-AIRPORT</td>\n",
       "      <td>core</td>\n",
       "      <td>CSD16</td>\n",
       "      <td>Lisbon</td>\n",
       "      <td>LIS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>CATCCGA</td>\n",
       "      <td>NaN</td>\n",
       "      <td>haib17KIU4866_H7HJMCCXY_SL272833</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SL272833</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>coastal</td>\n",
       "      <td>coastal</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17KIU4866_H7HJMCCXY_SL272830</th>\n",
       "      <td>CSD16-LIS-CTRL-LINEA</td>\n",
       "      <td>core</td>\n",
       "      <td>CSD16</td>\n",
       "      <td>Lisbon</td>\n",
       "      <td>LIS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>GTATACA</td>\n",
       "      <td>NaN</td>\n",
       "      <td>haib17KIU4866_H7HJMCCXY_SL272830</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SL272830</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>coastal</td>\n",
       "      <td>coastal</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17KIU4866_H7HJMCCXY_SL272831</th>\n",
       "      <td>CSD16-LIS-CTRL-LINEC</td>\n",
       "      <td>core</td>\n",
       "      <td>CSD16</td>\n",
       "      <td>Lisbon</td>\n",
       "      <td>LIS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>AGATCCC</td>\n",
       "      <td>NaN</td>\n",
       "      <td>haib17KIU4866_H7HJMCCXY_SL272831</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SL272831</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>coastal</td>\n",
       "      <td>coastal</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17KIU4866_H7HJMCCXY_SL272832</th>\n",
       "      <td>CSD16-LIS-CTRL-LINED</td>\n",
       "      <td>core</td>\n",
       "      <td>CSD16</td>\n",
       "      <td>Lisbon</td>\n",
       "      <td>LIS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>GCCGTAT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>haib17KIU4866_H7HJMCCXY_SL272832</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SL272832</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>coastal</td>\n",
       "      <td>coastal</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17KIU4866_HTWH2CCXY_SL364012</th>\n",
       "      <td>CSD16-MRS-0047</td>\n",
       "      <td>core</td>\n",
       "      <td>CSD16</td>\n",
       "      <td>Marseille</td>\n",
       "      <td>MRS</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>AATACGAA-AGGCCTGC</td>\n",
       "      <td>NaN</td>\n",
       "      <td>haib17KIU4866_HTWH2CCXY_SL364012</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>SL364012</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>coastal</td>\n",
       "      <td>coastal</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 42 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            metasub_name core_project project  \\\n",
       "uuid                                                                            \n",
       "haib17KIU4866_H7HJMCCXY_SL272833  CSD16-LIS-CTRL-AIRPORT         core   CSD16   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830    CSD16-LIS-CTRL-LINEA         core   CSD16   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831    CSD16-LIS-CTRL-LINEC         core   CSD16   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832    CSD16-LIS-CTRL-LINED         core   CSD16   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012          CSD16-MRS-0047         core   CSD16   \n",
       "\n",
       "                                       city city_code  latitude  longitude  \\\n",
       "uuid                                                                         \n",
       "haib17KIU4866_H7HJMCCXY_SL272833     Lisbon       LIS       NaN        NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830     Lisbon       LIS       NaN        NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831     Lisbon       LIS       NaN        NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832     Lisbon       LIS       NaN        NaN   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012  Marseille       MRS       NaN        NaN   \n",
       "\n",
       "                                 surface_material control_type elevation  ...  \\\n",
       "uuid                                                                      ...   \n",
       "haib17KIU4866_H7HJMCCXY_SL272833              NaN          NaN       NaN  ...   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830              NaN          NaN       NaN  ...   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831              NaN          NaN       NaN  ...   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832              NaN          NaN       NaN  ...   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012              NaN          NaN       NaN  ...   \n",
       "\n",
       "                                     index_sequence location_type  \\\n",
       "uuid                                                                \n",
       "haib17KIU4866_H7HJMCCXY_SL272833            CATCCGA           NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830            GTATACA           NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831            AGATCCC           NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832            GCCGTAT           NaN   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012  AATACGAA-AGGCCTGC           NaN   \n",
       "\n",
       "                                                  hudson_alpha_uid  \\\n",
       "uuid                                                                 \n",
       "haib17KIU4866_H7HJMCCXY_SL272833  haib17KIU4866_H7HJMCCXY_SL272833   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830  haib17KIU4866_H7HJMCCXY_SL272830   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831  haib17KIU4866_H7HJMCCXY_SL272831   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832  haib17KIU4866_H7HJMCCXY_SL272832   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012  haib17KIU4866_HTWH2CCXY_SL364012   \n",
       "\n",
       "                                  other_project_uid sample_type   sl_name  \\\n",
       "uuid                                                                        \n",
       "haib17KIU4866_H7HJMCCXY_SL272833                NaN         NaN  SL272833   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830                NaN         NaN  SL272830   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831                NaN         NaN  SL272831   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832                NaN         NaN  SL272832   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012                NaN         NaN  SL364012   \n",
       "\n",
       "                                  surface_ontology_fine  \\\n",
       "uuid                                                      \n",
       "haib17KIU4866_H7HJMCCXY_SL272833                    NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830                    NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831                    NaN   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832                    NaN   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012                    NaN   \n",
       "\n",
       "                                  surface_ontology_coarse  coastal  \\\n",
       "uuid                                                                 \n",
       "haib17KIU4866_H7HJMCCXY_SL272833                      NaN  coastal   \n",
       "haib17KIU4866_H7HJMCCXY_SL272830                      NaN  coastal   \n",
       "haib17KIU4866_H7HJMCCXY_SL272831                      NaN  coastal   \n",
       "haib17KIU4866_H7HJMCCXY_SL272832                      NaN  coastal   \n",
       "haib17KIU4866_HTWH2CCXY_SL364012                      NaN  coastal   \n",
       "\n",
       "                                  city_elevation  \n",
       "uuid                                              \n",
       "haib17KIU4866_H7HJMCCXY_SL272833         coastal  \n",
       "haib17KIU4866_H7HJMCCXY_SL272830         coastal  \n",
       "haib17KIU4866_H7HJMCCXY_SL272831         coastal  \n",
       "haib17KIU4866_H7HJMCCXY_SL272832         coastal  \n",
       "haib17KIU4866_HTWH2CCXY_SL364012         coastal  \n",
       "\n",
       "[5 rows x 42 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meta = data.meta\n",
    "taxa = data.wide_taxa.fillna(0)\n",
    "amrs = data.amr_genes\n",
    "meta.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>taxa_name</th>\n",
       "      <th>Acanthamoeba polyphaga mimivirus</th>\n",
       "      <th>Acaryochloris marina</th>\n",
       "      <th>Acetoanaerobium sticklandii</th>\n",
       "      <th>Acetobacter aceti</th>\n",
       "      <th>Acetobacter pasteurianus</th>\n",
       "      <th>Acetobacter persici</th>\n",
       "      <th>Acetobacter pomorum</th>\n",
       "      <th>Acetobacter senegalensis</th>\n",
       "      <th>Acetobacter sp. SLV-7</th>\n",
       "      <th>Acetobacter tropicalis</th>\n",
       "      <th>...</th>\n",
       "      <th>endosymbiont of Acanthamoeba sp. UWC8</th>\n",
       "      <th>endosymbiont of unidentified scaly snail isolate Monju</th>\n",
       "      <th>gamma proteobacterium HdN1</th>\n",
       "      <th>halophilic archaeon DL31</th>\n",
       "      <th>halophilic archaeon True-ADL</th>\n",
       "      <th>methanogenic archaeon ISO4-H5</th>\n",
       "      <th>secondary endosymbiont of Ctenarytaina eucalypti</th>\n",
       "      <th>secondary endosymbiont of Heteropsylla cubana</th>\n",
       "      <th>secondary endosymbiont of Trabutina mannipara</th>\n",
       "      <th>uncultured crAssphage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254769</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>28.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>52.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254770</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>38.0</td>\n",
       "      <td>133.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254771</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>339.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>28.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>87.0</td>\n",
       "      <td>440.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254772</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>69.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254773</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>28.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>29.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>118.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 4430 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "taxa_name                         Acanthamoeba polyphaga mimivirus  \\\n",
       "sample                                                               \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                               0.0   \n",
       "\n",
       "taxa_name                         Acaryochloris marina  \\\n",
       "sample                                                   \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                   0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                   0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                   0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                   0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                   0.0   \n",
       "\n",
       "taxa_name                         Acetoanaerobium sticklandii  \\\n",
       "sample                                                          \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                          0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                          0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                          0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                          0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                         28.0   \n",
       "\n",
       "taxa_name                         Acetobacter aceti  Acetobacter pasteurianus  \\\n",
       "sample                                                                          \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                5.0                      28.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770               38.0                     133.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771               48.0                     339.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772               24.0                      69.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773               30.0                      85.0   \n",
       "\n",
       "taxa_name                         Acetobacter persici  Acetobacter pomorum  \\\n",
       "sample                                                                       \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                  0.0                  0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                  0.0                  0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                 16.0                  0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                  0.0                  0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                 25.0                  3.0   \n",
       "\n",
       "taxa_name                         Acetobacter senegalensis  \\\n",
       "sample                                                       \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                       5.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                       5.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                      19.0   \n",
       "\n",
       "taxa_name                         Acetobacter sp. SLV-7  \\\n",
       "sample                                                    \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                    0.0   \n",
       "\n",
       "taxa_name                         Acetobacter tropicalis  ...  \\\n",
       "sample                                                    ...   \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                     0.0  ...   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                     0.0  ...   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                    28.0  ...   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                     0.0  ...   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                    29.0  ...   \n",
       "\n",
       "taxa_name                         endosymbiont of Acanthamoeba sp. UWC8  \\\n",
       "sample                                                                    \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                                    0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                                   14.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                                    0.0   \n",
       "\n",
       "taxa_name                         endosymbiont of unidentified scaly snail isolate Monju  \\\n",
       "sample                                                                                     \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                                               52.0        \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                                               37.0        \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                                               87.0        \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                                                0.0        \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                                               45.0        \n",
       "\n",
       "taxa_name                         gamma proteobacterium HdN1  \\\n",
       "sample                                                         \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                         0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                         0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                       440.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                        30.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                        23.0   \n",
       "\n",
       "taxa_name                         halophilic archaeon DL31  \\\n",
       "sample                                                       \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                       0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                       0.0   \n",
       "\n",
       "taxa_name                         halophilic archaeon True-ADL  \\\n",
       "sample                                                           \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                           0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                           0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                           0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                           0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                           0.0   \n",
       "\n",
       "taxa_name                         methanogenic archaeon ISO4-H5  \\\n",
       "sample                                                            \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                            0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                            0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                            0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                            0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                            0.0   \n",
       "\n",
       "taxa_name                         secondary endosymbiont of Ctenarytaina eucalypti  \\\n",
       "sample                                                                               \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                                               0.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                                               0.0   \n",
       "\n",
       "taxa_name                         secondary endosymbiont of Heteropsylla cubana  \\\n",
       "sample                                                                            \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                                           25.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                                           32.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                                           13.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                                            6.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                                            0.0   \n",
       "\n",
       "taxa_name                         secondary endosymbiont of Trabutina mannipara  \\\n",
       "sample                                                                            \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                                           14.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                                           21.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                                            8.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                                            6.0   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                                            0.0   \n",
       "\n",
       "taxa_name                         uncultured crAssphage  \n",
       "sample                                                   \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                    0.0  \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                    0.0  \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                    0.0  \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                    0.0  \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                  118.0  \n",
       "\n",
       "[5 rows x 4430 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "taxa.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4531, 1952)\n",
      "(4531, 3)\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>city_koppen_climate</th>\n",
       "      <th>setting</th>\n",
       "      <th>elevation</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254769</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254770</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254771</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254772</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254773</th>\n",
       "      <td>humid_subtropical_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                        city_koppen_climate  setting elevation\n",
       "sample                                                                        \n",
       "haib17CEM4890_H2NYMCCXY_SL254769   tropical_savanna_climate  unknown   unknown\n",
       "haib17CEM4890_H2NYMCCXY_SL254770   tropical_savanna_climate  unknown   unknown\n",
       "haib17CEM4890_H2NYMCCXY_SL254771   tropical_savanna_climate  unknown   unknown\n",
       "haib17CEM4890_H2NYMCCXY_SL254772   tropical_savanna_climate  unknown   unknown\n",
       "haib17CEM4890_H2NYMCCXY_SL254773  humid_subtropical_climate  unknown   unknown"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "design = meta.loc[taxa.index, ['city_koppen_climate', 'setting', 'elevation']].fillna('unknown')\n",
    "#design['city_population_density'] = design['city_population_density'].map(lambda x: float(x) if x != 'unknown' else float('nan'))\n",
    "#design['continent'] = [str(el) for el in design['continent']]\n",
    "#design['surface_ontology_fine'] = [str(el) for el in design['surface_ontology_fine']]\n",
    "#design = design.query('continent != \"0\"').query('city_population_density > 0')\n",
    "taxa = taxa.loc[design.index]\n",
    "\n",
    "taxa_prev = (taxa > 0).mean()\n",
    "taxa = taxa[taxa_prev[taxa_prev > 0.25].index]\n",
    "taxa += 1\n",
    "\n",
    "print(taxa.shape)\n",
    "print(design.shape)\n",
    "design.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading mavric\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Confirmed inputs\"\n",
      "[1] \"running HTS\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: converting counts to integer mode\n",
      "\n",
      "R[write to console]: -- note: fitType='parametric', but the dispersion trend was not well captured by the\n",
      "   function: y = a/x + b, and a local regression fit was automatically substituted.\n",
      "   specify fitType='local' or 'mean' to avoid this message next time.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Built PCA Obj\"\n",
      "[1] \"Finding relevant PCs\"\n"
     ]
    }
   ],
   "source": [
    "%%R -i taxa -i design\n",
    "\n",
    "devtools::load_all('/home/dcdanko/Dev/mavric/')\n",
    "\n",
    "mav = estVC(t(taxa), design, clist, verbose=T, autosel=F, sigcor=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                      setting\\n2% \n",
      "                                       0.00418853 \n",
      "                                    elevation\\n2% \n",
      "                                       0.03722007 \n",
      "              city_koppen_climate\\n6%&setting\\n2% \n",
      "                                       0.04271508 \n",
      "city_koppen_climate\\n6%&elevation\\n2%&setting\\n2% \n",
      "                                       0.55862373 \n",
      "                        elevation\\n2%&setting\\n2% \n",
      "                                       1.07463200 \n",
      "                          city_koppen_climate\\n6% \n",
      "                                       4.91026690 \n",
      "                                   Discarded\\n47% \n",
      "                                      46.56743136 \n",
      "                                 Unexplained\\n47% \n",
      "                                      46.80492232 \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "        \n",
    "design$city_koppen_climate = as.factor(design$city_koppen_climate)\n",
    "design$setting = as.factor(design$setting)\n",
    "design$elevation = as.factor(design$elevation)\n",
    "\n",
    "plotVars(mav, design, plotve=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['metasub_name',\n",
       " 'core_project',\n",
       " 'project',\n",
       " 'city',\n",
       " 'city_code',\n",
       " 'latitude',\n",
       " 'longitude',\n",
       " 'surface_material',\n",
       " 'control_type',\n",
       " 'elevation',\n",
       " 'line',\n",
       " 'station',\n",
       " 'surface',\n",
       " 'temperature',\n",
       " 'setting',\n",
       " 'num_reads',\n",
       " 'library_post_PCR_Qubit',\n",
       " 'library_QC_concentration',\n",
       " 'city_latitude',\n",
       " 'city_longitude',\n",
       " 'coastal_city',\n",
       " 'city_total_population',\n",
       " 'city_population_density',\n",
       " 'city_land_area_km2',\n",
       " 'city_ave_june_temp_c',\n",
       " 'city_elevation_meters',\n",
       " 'continent',\n",
       " 'city_koppen_climate',\n",
       " 'barcode',\n",
       " 'ha_id',\n",
       " 'hudson_alpha_flowcell',\n",
       " 'hudson_alpha_project',\n",
       " 'index_sequence',\n",
       " 'location_type',\n",
       " 'hudson_alpha_uid',\n",
       " 'other_project_uid',\n",
       " 'sample_type',\n",
       " 'sl_name',\n",
       " 'surface_ontology_fine',\n",
       " 'surface_ontology_coarse',\n",
       " 'coastal',\n",
       " 'city_elevation']"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(meta.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "amr_snames = amrs.index.to_list()\n",
    "amrs2 = (1000 * 1000 * amrs).applymap(int)\n",
    "amr_design =design.loc[amr_snames]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading mavric\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Confirmed inputs\"\n",
      "[1] \"running HTS\"\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc,  : \n",
      "  every gene contains at least one zero, cannot compute log geometric means\n",
      "Calls: <Anonymous> ... estimateSizeFactors -> .local -> estimateSizeFactorsForMatrix\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Error in estimateSizeFactorsForMatrix(counts(object), locfunc = locfunc,  : \n",
      "  every gene contains at least one zero, cannot compute log geometric means\n",
      "Calls: <Anonymous> ... estimateSizeFactors -> .local -> estimateSizeFactorsForMatrix\n"
     ]
    }
   ],
   "source": [
    "%%R -i amrs2 -i amr_design\n",
    "\n",
    "devtools::load_all('/home/dcdanko/Dev/mavric/')\n",
    "\n",
    "mav_amr = estVC(t(amrs2), amr_design, clist, verbose=T, autosel=F, sigcor=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "int() argument must be a string, a bytes-like object or a number, not 'DataFrame'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-fcdb8989a0ac>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mamrs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: int() argument must be a string, a bytes-like object or a number, not 'DataFrame'"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>city_koppen_climate</th>\n",
       "      <th>setting</th>\n",
       "      <th>elevation</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254769</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254770</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254771</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254772</th>\n",
       "      <td>tropical_savanna_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>haib17CEM4890_H2NYMCCXY_SL254773</th>\n",
       "      <td>humid_subtropical_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sossowski_BarcelonaNov2018_CSD16-BCN-258-29786-GGACTCCT-TATCCTCT</th>\n",
       "      <td>hot-summer_mediterranean_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sossowski_BarcelonaNov2018_CSD16-BCN-259-29787-GGACTCCT-GTAAGGAG</th>\n",
       "      <td>hot-summer_mediterranean_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sossowski_BarcelonaNov2018_CSD16-BCN-261-29787-GGACTCCT-AAGGAGTA</th>\n",
       "      <td>hot-summer_mediterranean_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sossowski_BarcelonaNov2018_CSD16-BCN-262-29787-GGACTCCT-CTAAGCCT</th>\n",
       "      <td>hot-summer_mediterranean_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sossowski_BarcelonaNov2018_CSD16-BCN-263-29787-GGACTCCT-CGTCTAAT</th>\n",
       "      <td>hot-summer_mediterranean_climate</td>\n",
       "      <td>unknown</td>\n",
       "      <td>unknown</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4531 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                                 city_koppen_climate  \\\n",
       "sample                                                                                 \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                            tropical_savanna_climate   \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                            tropical_savanna_climate   \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                            tropical_savanna_climate   \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                            tropical_savanna_climate   \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                           humid_subtropical_climate   \n",
       "...                                                                              ...   \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-258-29786-...  hot-summer_mediterranean_climate   \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-259-29787-...  hot-summer_mediterranean_climate   \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-261-29787-...  hot-summer_mediterranean_climate   \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-262-29787-...  hot-summer_mediterranean_climate   \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-263-29787-...  hot-summer_mediterranean_climate   \n",
       "\n",
       "                                                    setting elevation  \n",
       "sample                                                                 \n",
       "haib17CEM4890_H2NYMCCXY_SL254769                    unknown   unknown  \n",
       "haib17CEM4890_H2NYMCCXY_SL254770                    unknown   unknown  \n",
       "haib17CEM4890_H2NYMCCXY_SL254771                    unknown   unknown  \n",
       "haib17CEM4890_H2NYMCCXY_SL254772                    unknown   unknown  \n",
       "haib17CEM4890_H2NYMCCXY_SL254773                    unknown   unknown  \n",
       "...                                                     ...       ...  \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-258-29786-...  unknown   unknown  \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-259-29787-...  unknown   unknown  \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-261-29787-...  unknown   unknown  \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-262-29787-...  unknown   unknown  \n",
       "sossowski_BarcelonaNov2018_CSD16-BCN-263-29787-...  unknown   unknown  \n",
       "\n",
       "[4531 rows x 3 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [\n",
    "    'city',\n",
    "    'surface',\n",
    "    'surface_ontology_fine',\n",
    "    'city_koppen_climate',\n",
    "    'setting',\n",
    "    'elevation',\n",
    "    'coastal',\n",
    "    'hudson_alpha_flowcell',\n",
    "    'continent',\n",
    "]\n",
    "\n",
    "\n",
    "design = meta.loc[taxa.index, cols]\n",
    "\n",
    "for col in cols:\n",
    "    if col == 'continent':\n",
    "        continue\n",
    "    design[col] = design[col].fillna('unknown')\n",
    "\n",
    "def normalize_surface(el):\n",
    "    try:\n",
    "        el = el.lower()\n",
    "    except AttributeError:\n",
    "        return el\n",
    "    el = '_'.join(el.split())\n",
    "    return {\n",
    "         'ground': 'floor',\n",
    "         'palm_left': 'human_hand',\n",
    "         'vert_pole': 'pole',\n",
    "         'palm_right': 'human_hand',\n",
    "         'pedestrian_crossing_button': 'button',\n",
    "         'wooden bench': 'seat',\n",
    "         'gargabe_can': 'garbage',\n",
    "         'horiz_pole': 'railing',\n",
    "         'stairwell railing': 'railing',\n",
    "         'bike kiosk': 'kiosk',\n",
    "         'bench1': 'seat',\n",
    "         'kiosk2': 'kiosk',\n",
    "         'kiosk1': 'kiosk',\n",
    "         'ticket_machine': 'kiosk',\n",
    "         'escalator_handrail': 'railing',\n",
    "         'bench2': 'seat',\n",
    "         'lift_buttons': 'button',\n",
    "         'overhead_handrail': 'railing',\n",
    "         'poll': 'pole',\n",
    "        'garbage_can': 'garbage',\n",
    "        'bench': 'seat',\n",
    "        'door_handle;door_handle': 'door_knob',\n",
    "        'glass': 'window',\n",
    "        'ceiling_rail': 'railing',\n",
    "        'seat_rail': 'railing',\n",
    "        'ticket_hall;ticket_machine': 'kiosk',\n",
    "        'bench;platform': 'seat',\n",
    "        'platform;bench': 'seat',\n",
    "        'ticket_machine;ticket_hall': 'kiosk',\n",
    "        'wooden_bench': 'seat',\n",
    "        'stairwell_railing': 'railing',\n",
    "        'trolley_handle': '',\n",
    "        'turnstile_or_alternatives;turnstile_or_alternatives': 'turnstile',\n",
    "        'ticket_kiosks;ticket_kiosks': 'kiosk',\n",
    "        'ticket_machine;ticket_machine': 'kiosk',\n",
    "        'station_eletronic_kiosk': 'kiosk',\n",
    "        'station_railing': 'railing',\n",
    "        'station_seat': 'seat',\n",
    "        'center_seat': 'seat',\n",
    "        'negative_control_(air)': '',\n",
    "        'ceiling_rail;ceiling_rail': 'railing',\n",
    "        'rail': 'railing',\n",
    "        'seat_near_door': 'seat',\n",
    "        'ticketing_machine;ticketing_machine_': 'kiosk',\n",
    "        'jetway_1_seats_waiting_area': 'seat',\n",
    "    }.get(el, el)\n",
    "\n",
    "design['surface'] = design['surface'].map(normalize_surface)\n",
    "n_samps = design['surface'].value_counts()\n",
    "allowed_surfaces = set(n_samps[n_samps >= 8].index)\n",
    "design['surface'] = design['surface'].map(lambda el: 'other' if el not in allowed_surfaces else el)\n",
    "\n",
    "design = design.query('continent != \"Nan\"')\n",
    "design['continent'] = design['continent'].map(str)\n",
    "design = design.dropna()\n",
    "taxa = taxa.loc[design.index]\n",
    "\n",
    "taxa_prev = (taxa > 0).mean()\n",
    "taxa = taxa[taxa_prev[taxa_prev > 0.25].index]\n",
    "taxa += 1\n",
    "\n",
    "print(taxa.shape)\n",
    "print(design.shape)\n",
    "design.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "design['continent'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i taxa -i design\n",
    "\n",
    "devtools::load_all('/home/dcdanko/Dev/mavric/')\n",
    "\n",
    "mav = estVC(t(taxa), design, clist, verbose=T, autosel=F, sigcor=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i design\n",
    "\n",
    "attributes(mav)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "mav$pcavar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "plotVars <- function(results, annotation, plotzero = TRUE, incvals = FALSE, plotve = TRUE, maxStress = .1) {\n",
    "    results$pcs <- lapply(\n",
    "        results$pcs,\n",
    "        function(e) e[as.numeric(rownames(e)) > 0, , drop = F]\n",
    "    )\n",
    "    print(10)\n",
    "    vars <- round(\n",
    "        unlist(\n",
    "            sapply(1:length(results$pcs), function(p) {\n",
    "                if(nrow(results$pcs[[p]]) == 0) return(0)\n",
    "                return(\n",
    "                    sum(\n",
    "                        results$pcavar[as.numeric(rownames(results$pcs[[p]]))] *\n",
    "                        sapply(\n",
    "                            as.numeric(rownames(results$pcs[[p]])),\n",
    "                            function(i) totVarEx(results$pcasp, annotation, i, p)\n",
    "                        )\n",
    "                    )\n",
    "                )\n",
    "            })\n",
    "        )\n",
    "    )\n",
    "    print(20)\n",
    "    vn <- paste(colnames(annotation), paste(vars, \"%\", sep=''), sep='\\n')\n",
    "    print(vn)\n",
    "    print(21)\n",
    "    pcex <- sort(unique(unlist(lapply(results$pcs, rownames))))\n",
    "    print(22)\n",
    "    matrixed <- matrix(\n",
    "        unlist(\n",
    "            sapply(1:length(results$pcs), function(li) {\n",
    "                rv <- rep(0, length(pcex))\n",
    "                if(nrow(results$pcs[[li]]) > 0) {\n",
    "                    e <- rownames(results$pcs[[li]])\n",
    "                    for(i in 1:length(e)) rv[match(e[i], pcex)] <- results$pcavar[as.numeric(e[i])]*totVarEx(results$pcasp, annotation, as.numeric(e[i]), li)\n",
    "                }\n",
    "                return(rv)\n",
    "            })\n",
    "        ),\n",
    "        nrow=length(pcex),\n",
    "        dimnames=list(pcex, vn)\n",
    "    )\n",
    "    print(25)\n",
    "    applied <- apply(\n",
    "        matrixed,\n",
    "        1,\n",
    "        function(i) {\n",
    "            print(i)\n",
    "            print(25.1)\n",
    "            myorder <- order(i, decreasing=T)\n",
    "            print(25.15)\n",
    "            myn = sum((i > 0) & !is.na(i))\n",
    "            print(25.18)\n",
    "            print(myn)\n",
    "            print(myorder)\n",
    "            nz <- head(myorder, n=myn)\n",
    "            print(25.2)\n",
    "            if(length(nz) == 1) return(rbind(i[nz], vn[nz]))\n",
    "            print(25.3)\n",
    "            return(rbind(c(diff(t(embed(i[nz], 2))), i[nz[length(nz)]]), sapply(1:length(nz), function(n) paste(vn[nz[1:n]], collapse='&'))))\n",
    "        }\n",
    "    )\n",
    "    print(27)\n",
    "    unlisted <- unlist(applied)\n",
    "    print(28)\n",
    "    left_matrix <- matrix(unlisted, nrow=2)\n",
    "    print(29)\n",
    "    right_matrix <- matrix(\n",
    "        sapply(\n",
    "            unlist(\n",
    "                sapply(1:length(results$pcs), function(li) if(nrow(results$pcs[[li]]) == 0) return(c(0, vn[li])))\n",
    "            ),\n",
    "            function(i) i\n",
    "        ),\n",
    "        nrow=2\n",
    "    )\n",
    "    pcm <- cbind(left_matrix, right_matrix)\n",
    "    print(30)\n",
    "    unacc <- 100-sum(results$pcavar)\n",
    "    unex <- 100-(unacc+sum(as.numeric(pcm[1,])))\n",
    "    other.var <- c(unex, unacc)\n",
    "    names(other.var) <- paste(c(\"Unexplained\", \"Discarded\"), paste(as.character(round(c(unex, unacc))), \"%\", sep=''), sep='\\n')\n",
    "    pl.vars <- sort(c(tapply(as.numeric(unlist(pcm[1,])), unlist(lapply(lapply(strsplit(unlist(pcm[2,]), \"&\"), sort), paste, collapse = '&')), sum), other.var), decreasing = incvals)\n",
    "    if(!plotzero) pl.vars <- pl.vars[pl.vars > 0]\n",
    "    if(!plotve) return(pl.vars)\n",
    "    pl.vars.c <- ceiling(pl.vars)\n",
    "    repeat {\n",
    "        pl.tmp <- abs(jitter(pl.vars.c))\n",
    "        if(plotzero) pl.tmp[pl.vars.c == 0] = 0\n",
    "        pl.ve <- venneuler::venneuler(pl.tmp)\n",
    "        if(pl.ve$stress < maxStress) break\n",
    "    }\n",
    "    print(40)\n",
    "    if(plotve) plot(venneuler::venneuler(pl.tmp))\n",
    "    return(pl.vars)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i design\n",
    "\n",
    "design$city_koppen_climate = as.factor(design$city_koppen_climate)\n",
    "design$setting = as.factor(design$setting)\n",
    "design$elevation = as.factor(design$elevation)\n",
    "design$city = as.factor(design$city)\n",
    "design$surface = as.factor(design$surface)\n",
    "design$surface_ontology_fine = as.factor(design$surface_ontology_fine)\n",
    "design$coastal = as.factor(design$coastal)\n",
    "design$hudson_alpha_flowcell = as.factor(design$hudson_alpha_flowcell)\n",
    "#design$continent = as.factor(design$continent)\n",
    "\n",
    "plotVars(mav, design, plotve=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [\n",
    "    'city',\n",
    "    'surface',\n",
    "    'surface_ontology_fine',\n",
    "    'city_koppen_climate',\n",
    "    'setting',\n",
    "    'elevation',\n",
    "    'coastal',\n",
    "    'hudson_alpha_flowcell',\n",
    "]\n",
    "\n",
    "\n",
    "design = meta.loc[taxa.index, cols]\n",
    "\n",
    "for col in cols:\n",
    "    if col == 'continent':\n",
    "        continue\n",
    "    design[col] = design[col].fillna('unknown')\n",
    "\n",
    "def normalize_surface(el):\n",
    "    try:\n",
    "        el = el.lower()\n",
    "    except AttributeError:\n",
    "        return el\n",
    "    el = '_'.join(el.split())\n",
    "    return {\n",
    "         'ground': 'floor',\n",
    "         'palm_left': 'human_hand',\n",
    "         'vert_pole': 'pole',\n",
    "         'palm_right': 'human_hand',\n",
    "         'pedestrian_crossing_button': 'button',\n",
    "         'wooden bench': 'seat',\n",
    "         'gargabe_can': 'garbage',\n",
    "         'horiz_pole': 'railing',\n",
    "         'stairwell railing': 'railing',\n",
    "         'bike kiosk': 'kiosk',\n",
    "         'bench1': 'seat',\n",
    "         'kiosk2': 'kiosk',\n",
    "         'kiosk1': 'kiosk',\n",
    "         'ticket_machine': 'kiosk',\n",
    "         'escalator_handrail': 'railing',\n",
    "         'bench2': 'seat',\n",
    "         'lift_buttons': 'button',\n",
    "         'overhead_handrail': 'railing',\n",
    "         'poll': 'pole',\n",
    "        'garbage_can': 'garbage',\n",
    "        'bench': 'seat',\n",
    "        'door_handle;door_handle': 'door_knob',\n",
    "        'glass': 'window',\n",
    "        'ceiling_rail': 'railing',\n",
    "        'seat_rail': 'railing',\n",
    "        'ticket_hall;ticket_machine': 'kiosk',\n",
    "        'bench;platform': 'seat',\n",
    "        'platform;bench': 'seat',\n",
    "        'ticket_machine;ticket_hall': 'kiosk',\n",
    "        'wooden_bench': 'seat',\n",
    "        'stairwell_railing': 'railing',\n",
    "        'trolley_handle': '',\n",
    "        'turnstile_or_alternatives;turnstile_or_alternatives': 'turnstile',\n",
    "        'ticket_kiosks;ticket_kiosks': 'kiosk',\n",
    "        'ticket_machine;ticket_machine': 'kiosk',\n",
    "        'station_eletronic_kiosk': 'kiosk',\n",
    "        'station_railing': 'railing',\n",
    "        'station_seat': 'seat',\n",
    "        'center_seat': 'seat',\n",
    "        'negative_control_(air)': '',\n",
    "        'ceiling_rail;ceiling_rail': 'railing',\n",
    "        'rail': 'railing',\n",
    "        'seat_near_door': 'seat',\n",
    "        'ticketing_machine;ticketing_machine_': 'kiosk',\n",
    "        'jetway_1_seats_waiting_area': 'seat',\n",
    "    }.get(el, el)\n",
    "\n",
    "design['surface'] = design['surface'].map(normalize_surface)\n",
    "n_samps = design['surface'].value_counts()\n",
    "allowed_surfaces = set(n_samps[n_samps >= 8].index)\n",
    "design['surface'] = design['surface'].map(lambda el: 'other' if el not in allowed_surfaces else el)\n",
    "\n",
    "\n",
    "\n",
    "taxa = taxa.loc[design.index]\n",
    "\n",
    "taxa_prev = (taxa > 0).mean()\n",
    "taxa = taxa[taxa_prev[taxa_prev > 0.25].index]\n",
    "taxa += 1\n",
    "\n",
    "print(taxa.shape)\n",
    "print(design.shape)\n",
    "design.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i taxa -i design\n",
    "\n",
    "devtools::load_all('/home/dcdanko/Dev/mavric/')\n",
    "\n",
    "mav = estVC(t(taxa), design, clist, verbose=T, autosel=F, sigcor=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "design$city_koppen_climate = as.factor(design$city_koppen_climate)\n",
    "design$setting = as.factor(design$setting)\n",
    "design$elevation = as.factor(design$elevation)\n",
    "design$city = as.factor(design$city)\n",
    "design$surface = as.factor(design$surface)\n",
    "design$surface_ontology_fine = as.factor(design$surface_ontology_fine)\n",
    "design$coastal = as.factor(design$coastal)\n",
    "design$hudson_alpha_flowcell = as.factor(design$hudson_alpha_flowcell)\n",
    "\n",
    "plotVars(mav, design, plotve=F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}