8.gtex-interpret/0.process-gtex-blood-vae-feature.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Capturing Blood Tissue in GTEx Data\n",
    "\n",
    "**Gregory Way, 2018**\n",
    "\n",
    "When viewing sample correlation differences across z dimensions stratified by tissue-type, we often observed a rapid increase in correlation after increasing model capacity by one.\n",
    "\n",
    "For example, the ability of variational autoencoders to capture blood tissue correlation jumps by nearly 50% between bottleneck dimensions 2 and 3 (see below).\n",
    "\n",
    "![sample-correlation_Blood_GTEX_signal_pearson.png](https://raw.githubusercontent.com/greenelab/BioBombe/master/4.analyze-components/figures/GTEX/sample-correlation/sample-type/sample-correlation_Blood_GTEX_signal_pearson.png)\n",
    "\n",
    "## Procedure\n",
    "\n",
    "In the following notebook, we extract two representative weight matrices for VAE latent space dimensions 2 and 3.\n",
    "We apply two compression feature interpretaion approaches to the weight vectors for each latent space feature. \n",
    "\n",
    "1. Our matrix interpretation approach\n",
    "2. Overrepresentation Tests using high weight genes\n",
    "\n",
    "In both approaches we use genesets derived in the xCell paper that represent cell-types ([Aran et al. 2017](https://doi.org/10.1186/s13059-017-1349-1))\n",
    "\n",
    "We output the results of both approaches and analyze the results in subsequent notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import glob\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "sys.path.append('../scripts')\n",
    "from latent import latentModel, load_hetnets, parse_gmt, run_overrepresentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "metaedge = 'GpXCELL'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_dir = os.path.join('..', '2.ensemble-z-analysis', 'results',\n",
    "                        'GTEX_results', 'ensemble_z_matrices')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_weight_2 = os.path.join(base_dir, 'gtex_components_2', 'model_908341_weight_matrix.tsv.gz')\n",
    "sample_weight_3 = os.path.join(base_dir, 'gtex_components_3', 'model_908341_weight_matrix.tsv.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>vae_0_two</th>\n",
       "      <th>vae_1_two</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3043</th>\n",
       "      <td>0.176201</td>\n",
       "      <td>-0.165548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3040</th>\n",
       "      <td>0.163633</td>\n",
       "      <td>-0.163657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6280</th>\n",
       "      <td>-0.007069</td>\n",
       "      <td>-0.225171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5644</th>\n",
       "      <td>0.219859</td>\n",
       "      <td>-0.106132</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2512</th>\n",
       "      <td>-0.059835</td>\n",
       "      <td>0.051747</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         vae_0_two  vae_1_two\n",
       "gene_id                      \n",
       "3043      0.176201  -0.165548\n",
       "3040      0.163633  -0.163657\n",
       "6280     -0.007069  -0.225171\n",
       "5644      0.219859  -0.106132\n",
       "2512     -0.059835   0.051747"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_weight_2_df = pd.read_table(sample_weight_2, index_col=0)\n",
    "sample_weight_3_df = pd.read_table(sample_weight_3, index_col=0)\n",
    "\n",
    "sample_weight_2_df = sample_weight_2_df.loc[:, sample_weight_2_df.columns.str.contains('vae')]\n",
    "sample_weight_3_df = sample_weight_3_df.loc[:, sample_weight_3_df.columns.str.contains('vae')]\n",
    "\n",
    "# Recode column names\n",
    "sample_weight_2_df.columns = sample_weight_2_df.columns + '_two'\n",
    "sample_weight_3_df.columns = sample_weight_3_df.columns + '_three'\n",
    "\n",
    "sample_weight_2_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>vae_0_two</th>\n",
       "      <th>vae_1_two</th>\n",
       "      <th>vae_0_three</th>\n",
       "      <th>vae_1_three</th>\n",
       "      <th>vae_2_three</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_id</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>3043</th>\n",
       "      <td>0.176201</td>\n",
       "      <td>-0.165548</td>\n",
       "      <td>0.314363</td>\n",
       "      <td>-0.170519</td>\n",
       "      <td>0.021420</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3040</th>\n",
       "      <td>0.163633</td>\n",
       "      <td>-0.163657</td>\n",
       "      <td>0.294492</td>\n",
       "      <td>-0.188674</td>\n",
       "      <td>0.008922</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         vae_0_two  vae_1_two  vae_0_three  vae_1_three  vae_2_three\n",
       "gene_id                                                             \n",
       "3043      0.176201  -0.165548     0.314363    -0.170519     0.021420\n",
       "3040      0.163633  -0.163657     0.294492    -0.188674     0.008922"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "combined_model = sample_weight_2_df.merge(sample_weight_3_df,\n",
    "                                          left_index=True,\n",
    "                                          right_index=True)\n",
    "combined_model.index = combined_model.index.astype('str')\n",
    "combined_model.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load hetnets for the given metaedge for rapid latent feature interpretation\n",
    "hetnets = load_hetnets(\n",
    "    hetnet_file='../3.build-hetnets/hetnets/interpret_hetnet.json.bz2',\n",
    "    permuted_directory='../3.build-hetnets/hetnets/permuted/',\n",
    "    subset_genes=combined_model.index,\n",
    "    metaedge_abbrev=metaedge\n",
    "   )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply Interpret Compression Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.3 s, sys: 18.2 ms, total: 6.32 s\n",
      "Wall time: 8.86 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "mult_results = {}\n",
    "all_list = []\n",
    "for model in hetnets.keys():\n",
    "    hetnet = hetnets[model]\n",
    "    mult_results[model] = combined_model.T @ hetnet\n",
    "    long_result = mult_results[model].reset_index().melt(id_vars=['index'])\n",
    "    long_result = long_result.assign(model=model)\n",
    "    if model != 'real':\n",
    "        long_result = long_result.assign(model_type='permuted')\n",
    "    else:\n",
    "        long_result = long_result.assign(model_type='real')\n",
    "    all_list.append(long_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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>full_feature</th>\n",
       "      <th>variable</th>\n",
       "      <th>raw_score</th>\n",
       "      <th>model</th>\n",
       "      <th>model_type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>-0.580380</td>\n",
       "      <td>0</td>\n",
       "      <td>permuted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>vae_1_two</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>-0.092745</td>\n",
       "      <td>0</td>\n",
       "      <td>permuted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>-0.390124</td>\n",
       "      <td>0</td>\n",
       "      <td>permuted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>vae_1_three</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>-0.430026</td>\n",
       "      <td>0</td>\n",
       "      <td>permuted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>vae_2_three</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>-0.081764</td>\n",
       "      <td>0</td>\n",
       "      <td>permuted</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  full_feature             variable  raw_score model model_type\n",
       "0    vae_0_two  Adipocytes_ENCODE_1  -0.580380     0   permuted\n",
       "1    vae_1_two  Adipocytes_ENCODE_1  -0.092745     0   permuted\n",
       "2  vae_0_three  Adipocytes_ENCODE_1  -0.390124     0   permuted\n",
       "3  vae_1_three  Adipocytes_ENCODE_1  -0.430026     0   permuted\n",
       "4  vae_2_three  Adipocytes_ENCODE_1  -0.081764     0   permuted"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_df = pd.concat(all_list)\n",
    "all_df.value = all_df.value.astype(float)\n",
    "all_df = all_df.rename(columns={'index': 'full_feature', 'value': 'raw_score'})\n",
    "all_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_group = all_df.groupby(['model_type', 'full_feature', 'variable'])\n",
    "all_group_mean = all_group.mean().reset_index()\n",
    "\n",
    "real_df = all_group_mean.query(\"model_type == 'real'\")\n",
    "permuted_mean = all_group_mean.query(\"model_type == 'permuted'\")\n",
    "\n",
    "permuted_std = (\n",
    "    all_group\n",
    "    .std()\n",
    "    .reset_index()\n",
    "    .query(\"model_type == 'permuted'\")\n",
    ")"
   ]
  },
  {
   "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></th>\n",
       "      <th>model_type</th>\n",
       "      <th>full_feature</th>\n",
       "      <th>variable</th>\n",
       "      <th>raw_score</th>\n",
       "      <th>z_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>real</td>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_ENCODE_1</td>\n",
       "      <td>0.251069</td>\n",
       "      <td>2.423840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>real</td>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_ENCODE_2</td>\n",
       "      <td>0.265078</td>\n",
       "      <td>2.788173</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>real</td>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_ENCODE_3</td>\n",
       "      <td>0.314689</td>\n",
       "      <td>1.930844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>real</td>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_FANTOM_1</td>\n",
       "      <td>0.437572</td>\n",
       "      <td>3.728542</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>real</td>\n",
       "      <td>vae_0_three</td>\n",
       "      <td>Adipocytes_FANTOM_2</td>\n",
       "      <td>0.437572</td>\n",
       "      <td>2.318366</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  model_type full_feature             variable  raw_score   z_score\n",
       "0       real  vae_0_three  Adipocytes_ENCODE_1   0.251069  2.423840\n",
       "1       real  vae_0_three  Adipocytes_ENCODE_2   0.265078  2.788173\n",
       "2       real  vae_0_three  Adipocytes_ENCODE_3   0.314689  1.930844\n",
       "3       real  vae_0_three  Adipocytes_FANTOM_1   0.437572  3.728542\n",
       "4       real  vae_0_three  Adipocytes_FANTOM_2   0.437572  2.318366"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z_score = (real_df.reset_index(drop=True).raw_score - permuted_mean.raw_score) / permuted_std.raw_score\n",
    "real_df = real_df.reset_index(drop=True).assign(z_score=z_score)\n",
    "\n",
    "real_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output z scores per xCell genesets for all VAE features (5 total)\n",
    "file = os.path.join('results', 'gtex_vae_example_interpret_compression.tsv')\n",
    "real_df.sort_values(by='z_score').to_csv(file, sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine the most distinguishing features between the two models. This will test\n",
    "# which genesets are the most differently enriched _in sum_ between z = 2 and z = 3\n",
    "feature_info_df = (\n",
    "    pd.DataFrame(real_df['full_feature'].str.split('_').values.tolist(),\n",
    "                 columns=['algorithm', 'feature', 'model_z'])\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_info_df = pd.concat([real_df, feature_info_df], axis=1)\n",
    "feature_info_df = feature_info_df.assign(abs_z_score = feature_info_df.z_score.abs())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>model_z</th>\n",
       "      <th>three</th>\n",
       "      <th>two</th>\n",
       "      <th>abs_diff</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>variable</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Skeletal-muscle_FANTOM_2</th>\n",
       "      <td>7.946984</td>\n",
       "      <td>15.068043</td>\n",
       "      <td>7.121059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Monocytes_FANTOM_2</th>\n",
       "      <td>6.910606</td>\n",
       "      <td>0.709876</td>\n",
       "      <td>6.200731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Monocytes_IRIS_3</th>\n",
       "      <td>7.092820</td>\n",
       "      <td>1.342540</td>\n",
       "      <td>5.750279</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Astrocytes_FANTOM_3</th>\n",
       "      <td>4.243827</td>\n",
       "      <td>9.696995</td>\n",
       "      <td>5.453168</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Monocytes_NOVERSHTERN_1</th>\n",
       "      <td>7.413432</td>\n",
       "      <td>1.991430</td>\n",
       "      <td>5.422002</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "model_z                      three        two  abs_diff\n",
       "variable                                               \n",
       "Skeletal-muscle_FANTOM_2  7.946984  15.068043  7.121059\n",
       "Monocytes_FANTOM_2        6.910606   0.709876  6.200731\n",
       "Monocytes_IRIS_3          7.092820   1.342540  5.750279\n",
       "Astrocytes_FANTOM_3       4.243827   9.696995  5.453168\n",
       "Monocytes_NOVERSHTERN_1   7.413432   1.991430  5.422002"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "feature_info_df = (\n",
    "    feature_info_df.groupby(['variable', 'model_z'])\n",
    "    .mean()\n",
    "    .reset_index()\n",
    "    .pivot(index='variable', columns='model_z', values = 'abs_z_score')\n",
    ")\n",
    "\n",
    "feature_info_df = (\n",
    "    feature_info_df.assign(abs_diff = (feature_info_df.three - feature_info_df.two).abs())\n",
    "    .sort_values(by='abs_diff', ascending=False)\n",
    ")\n",
    "\n",
    "feature_info_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output the most distinguishing features considering all features within each z dimension\n",
    "file = os.path.join('results', 'gtex_vae_example_differentiating_features.tsv')\n",
    "feature_info_df.to_csv(file, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform overrepresentation tests on the same features\n",
    "\n",
    "### Split into positive and negative tails and extract high weight genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "489"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geneset_file = os.path.join('..', '3.build-hetnets', 'data', 'xcell_all_entrez.gmt')\n",
    "xcell_genesets_gmt = parse_gmt(gene_sets=[geneset_file])\n",
    "len(xcell_genesets_gmt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "vae_2_lm = latentModel(filename=sample_weight_2,\n",
    "                       z_dim=2,\n",
    "                       dataset_name='GTEX',\n",
    "                       algorithm_name='VAE',\n",
    "                       weight_seed='908341',\n",
    "                       shuffled_true=False)\n",
    "\n",
    "vae_3_lm = latentModel(filename=sample_weight_3,\n",
    "                       z_dim=3,\n",
    "                       dataset_name='GTEX',\n",
    "                       algorithm_name='VAE',\n",
    "                       weight_seed='908341',\n",
    "                       shuffled_true=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "vae_2_lm.get_high_weight_genes()\n",
    "vae_3_lm.get_high_weight_genes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "background_genes = []\n",
    "for xcell_name, xcell_gene_set in xcell_genesets_gmt.items():\n",
    "    background_genes += xcell_gene_set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.56 s, sys: 18.4 ms, total: 4.57 s\n",
      "Wall time: 6.79 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "gene_list = vae_2_lm.w_df[vae_2_lm.pos_high_w_df.vae_0].index.tolist()\n",
    "vae_feature_0_zdim_2_pos = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_2_lm.w_df[vae_2_lm.neg_high_w_df.vae_0].index.tolist()\n",
    "vae_feature_0_zdim_2_neg = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_2_lm.w_df[vae_2_lm.pos_high_w_df.vae_1].index.tolist()\n",
    "vae_feature_1_zdim_2_pos = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_2_lm.w_df[vae_2_lm.neg_high_w_df.vae_1].index.tolist()\n",
    "vae_feature_1_zdim_2_neg = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.pos_high_w_df.vae_0].index.tolist()\n",
    "vae_feature_0_zdim_3_pos = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.neg_high_w_df.vae_0].index.tolist()\n",
    "vae_feature_0_zdim_3_neg = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.pos_high_w_df.vae_1].index.tolist()\n",
    "vae_feature_1_zdim_3_pos = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.neg_high_w_df.vae_1].index.tolist()\n",
    "vae_feature_1_zdim_3_neg = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.pos_high_w_df.vae_2].index.tolist()\n",
    "vae_feature_2_zdim_3_pos = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)\n",
    "\n",
    "gene_list = vae_3_lm.w_df[vae_3_lm.neg_high_w_df.vae_2].index.tolist()\n",
    "vae_feature_2_zdim_3_neg = run_overrepresentation(gene_list=gene_list,\n",
    "                                                  gene_set_dict=xcell_genesets_gmt,\n",
    "                                                  background_genes=background_genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "vae_feature_0_zdim_2_pos = vae_feature_0_zdim_2_pos.assign(feature='vae_0_two', tailed='pos')\n",
    "vae_feature_0_zdim_2_neg = vae_feature_0_zdim_2_neg.assign(feature='vae_0_two', tailed='neg')\n",
    "vae_feature_1_zdim_2_pos = vae_feature_1_zdim_2_pos.assign(feature='vae_1_two', tailed='pos')\n",
    "vae_feature_1_zdim_2_neg = vae_feature_1_zdim_2_neg.assign(feature='vae_1_two', tailed='neg')\n",
    "\n",
    "vae_feature_0_zdim_3_pos = vae_feature_0_zdim_3_pos.assign(feature='vae_0_three', tailed='pos')\n",
    "vae_feature_0_zdim_3_neg = vae_feature_0_zdim_3_neg.assign(feature='vae_0_three', tailed='neg')\n",
    "vae_feature_1_zdim_3_pos = vae_feature_1_zdim_3_pos.assign(feature='vae_1_three', tailed='pos')\n",
    "vae_feature_1_zdim_3_neg = vae_feature_1_zdim_3_neg.assign(feature='vae_1_three', tailed='neg')\n",
    "vae_feature_2_zdim_3_pos = vae_feature_2_zdim_3_pos.assign(feature='vae_2_three', tailed='pos')\n",
    "vae_feature_2_zdim_3_neg = vae_feature_2_zdim_3_neg.assign(feature='vae_2_three', tailed='neg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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>pval</th>\n",
       "      <th>odds</th>\n",
       "      <th>feature</th>\n",
       "      <th>tailed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>pro-B-cells_NOVERSHTERN_2</th>\n",
       "      <td>0.220006</td>\n",
       "      <td>4.339523</td>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>pos</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Endothelial-cells_ENCODE_3</th>\n",
       "      <td>0.369979</td>\n",
       "      <td>2.251034</td>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>pos</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CD4+-T-cells_FANTOM_3</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.732759</td>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>pos</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>DC_FANTOM_3</th>\n",
       "      <td>0.323462</td>\n",
       "      <td>2.681992</td>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>pos</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MPP_BLUEPRINT_1</th>\n",
       "      <td>0.199574</td>\n",
       "      <td>2.519897</td>\n",
       "      <td>vae_0_two</td>\n",
       "      <td>pos</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                pval      odds    feature tailed\n",
       "pro-B-cells_NOVERSHTERN_2   0.220006  4.339523  vae_0_two    pos\n",
       "Endothelial-cells_ENCODE_3  0.369979  2.251034  vae_0_two    pos\n",
       "CD4+-T-cells_FANTOM_3       1.000000  0.732759  vae_0_two    pos\n",
       "DC_FANTOM_3                 0.323462  2.681992  vae_0_two    pos\n",
       "MPP_BLUEPRINT_1             0.199574  2.519897  vae_0_two    pos"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overrepresented_results_df = pd.concat([\n",
    "    vae_feature_0_zdim_2_pos,\n",
    "    vae_feature_0_zdim_2_neg,\n",
    "    vae_feature_1_zdim_2_pos,\n",
    "    vae_feature_1_zdim_2_neg,\n",
    "    vae_feature_0_zdim_3_pos,\n",
    "    vae_feature_0_zdim_3_neg,\n",
    "    vae_feature_1_zdim_3_pos,\n",
    "    vae_feature_1_zdim_3_neg,\n",
    "    vae_feature_2_zdim_3_pos,\n",
    "    vae_feature_2_zdim_3_neg\n",
    "])\n",
    "\n",
    "overrepresented_results_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = os.path.join('results', 'gtex_vae_example_overrepresentation.tsv')\n",
    "(\n",
    "    overrepresented_results_df\n",
    "    .reset_index()\n",
    "    .rename(columns={'index': 'variable',\n",
    "                     'feature': 'full_feature'})\n",
    "    .sort_values(by='pval')\n",
    "    .to_csv(file, sep='\\t', index=False)\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:biobombe]",
   "language": "python",
   "name": "conda-env-biobombe-py"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}