genome_microbiome_interactions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "urban-windsor",
   "metadata": {},
   "source": [
    "# Genome-microbiome interactions\n",
    "\n",
    "We now will look at associations between microbes and metabolites that ar emodulted by the genetic background of the host. This essentially comes down to looking at interaction terms between genetics and the microbiome in a linear model with the metabolite as the response. We will again use the confounder-corrected metabolite abundances and do a train and validation split.\n",
    "\n",
    "Let's start by loading the required packages and reading the required data sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "canadian-camera",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyplink import PyPlink\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats.multitest import fdrcorrection\n",
    "from rich.progress import track\n",
    "import arivale_data_interface as adi\n",
    "from utils import rsid2gene\n",
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")\n",
    "\n",
    "metabolites = pd.read_csv(\"data/metabolites_residuals.csv\")\n",
    "microbes = pd.read_csv(\"data/genera_clr_filtered.csv\")\n",
    "microbes = microbes[microbes.stool_sample_id.isin(metabolites.stool_sample_id)]\n",
    "genotypes = PyPlink(\"input_bed/all_chr/all_genomes_09112019_all_chr\")\n",
    "\n",
    "# Add gene annotations to the variants by querying dbSNP\n",
    "sig_snps = pd.read_csv(\"data/final_results.csv\")\n",
    "sig_snps[\"rsid\"] = sig_snps.SNP.str.split(\";\").str[0]\n",
    "genes = rsid2gene(sig_snps.rsid.unique())\n",
    "sig_snps = sig_snps.merge(genes, left_on=\"rsid\", right_on=\"snp\", how=\"left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "leading-south",
   "metadata": {},
   "source": [
    "Now we can start with the interaction terms. We will run tests for all metabolites that have a significant genetic association. We start by assembling the metabolome and metabolite data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "revised-nashville",
   "metadata": {},
   "outputs": [],
   "source": [
    "metabolites_and_microbes = pd.merge(metabolites, microbes, on=\"stool_sample_id\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "downtown-updating",
   "metadata": {},
   "source": [
    "And we can run the association tests. For that we start by writing a function that runs the test for a single metabolite-microbe combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "controlling-guest",
   "metadata": {},
   "outputs": [],
   "source": [
    "MIN_N = 30\n",
    "genome_ids = genotypes.get_fam().iid.values\n",
    "\n",
    "def interactions(args):\n",
    "    \"\"\"Run the interaction analysis for a metabolite-microbe pair.\"\"\"\n",
    "    met, mic = args\n",
    "    snps = sig_snps.SNP[sig_snps.metabolite == met]\n",
    "    geno = pd.DataFrame({s: genotypes.get_geno_marker(s) for s in snps}, index=genome_ids)\n",
    "    df = pd.merge(metabolites[[\"genome_id\", \"stool_sample_id\", met]], geno, left_on=\"genome_id\", right_index=True)\n",
    "    df = pd.merge(df, metabolites_and_microbes[[\"stool_sample_id\"] + [mic]], on=\"stool_sample_id\").dropna(subset=[mic])\n",
    "    \n",
    "    result = pd.DataFrame(\n",
    "        index = [met + ':' + mic + ':' + i for i in snps],\n",
    "        columns = ['metabolite','taxon','snps', 'snp_beta', 'interaction_beta','p', \n",
    "                   'baseline_r2', 'interaction_r2', 'full_r2', 'n_major', 'n_heterozygous', 'n_homozygous']\n",
    "    )\n",
    "    \n",
    "    if df.shape[0] < MIN_N:\n",
    "        return result\n",
    "    \n",
    "    for snp in snps:\n",
    "        n_allele = df[snp].value_counts()\n",
    "        if (n_allele < MIN_N).any() or len(n_allele) < 3:\n",
    "            continue\n",
    "        gxe_term = f\"Q('{snp}'):Q('{mic}')\"\n",
    "        formula_base = f\"{met} ~ Q('{snp}') + Q('{mic}')\"\n",
    "        formula_full = f\"{met} ~ Q('{snp}') + Q('{mic}') + {gxe_term}\"\n",
    "        base_fit = smf.ols(formula_base, data=df).fit()\n",
    "        full_fit = smf.ols(formula_full, data=df).fit()\n",
    "        result.loc[met + ':' + mic + ':' + snp] = [\n",
    "            met, \n",
    "            mic, \n",
    "            snp,\n",
    "            full_fit.params.loc[f\"Q('{snp}')\"],\n",
    "            full_fit.params.loc[gxe_term],\n",
    "            full_fit.pvalues.loc[gxe_term],\n",
    "            base_fit.rsquared,\n",
    "            full_fit.rsquared - base_fit.rsquared,\n",
    "            full_fit.rsquared,\n",
    "            n_allele.get(0, 0),\n",
    "            n_allele.get(1, 0),\n",
    "            n_allele.get(2, 0)\n",
    "        ]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "trained-enlargement",
   "metadata": {},
   "source": [
    "With that we can run our models for all combinations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "isolated-martin",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2cd62d3aab284e54b92c657a78dbfe4c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from os import path\n",
    "from multiprocessing import Pool\n",
    "\n",
    "if not path.exists(\"data/interaction_results1.csv\"):\n",
    "    with Pool(8) as pool:\n",
    "        args = [(met, mic) for met in list(sig_snps.metabolite.unique()) for mic in microbes.columns[1:]]\n",
    "        it = track(pool.imap_unordered(interactions, args), total=len(args), description='Fitting models')\n",
    "        results = pd.concat(list(it))\n",
    "    results['q'] = fdrcorrection(results.p)[1]\n",
    "    results.to_csv(\"data/interaction_results.csv\", index=False)\n",
    "else:\n",
    "    results = pd.read_csv(\"data/interaction_results.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bronze-assumption",
   "metadata": {},
   "source": [
    "Now we merge in the annotations for genes and metabolites and extract the bacterial genus and family."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "concerned-consequence",
   "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>metabolite</th>\n",
       "      <th>taxon</th>\n",
       "      <th>snps</th>\n",
       "      <th>snp_beta</th>\n",
       "      <th>interaction_beta</th>\n",
       "      <th>p</th>\n",
       "      <th>baseline_r2</th>\n",
       "      <th>interaction_r2</th>\n",
       "      <th>full_r2</th>\n",
       "      <th>n_major</th>\n",
       "      <th>...</th>\n",
       "      <th>SUPER_PATHWAY</th>\n",
       "      <th>SUB_PATHWAY</th>\n",
       "      <th>CAS</th>\n",
       "      <th>KEGG</th>\n",
       "      <th>SNP</th>\n",
       "      <th>CHR</th>\n",
       "      <th>genes</th>\n",
       "      <th>rsid</th>\n",
       "      <th>genus</th>\n",
       "      <th>family</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9687</th>\n",
       "      <td>metabolite_100001104</td>\n",
       "      <td>Ruminococcaceae|Ruminococcaceae_UCG-002</td>\n",
       "      <td>rs56672945</td>\n",
       "      <td>0.309384</td>\n",
       "      <td>-0.032553</td>\n",
       "      <td>0.000001</td>\n",
       "      <td>0.189084</td>\n",
       "      <td>0.014297</td>\n",
       "      <td>0.203381</td>\n",
       "      <td>944</td>\n",
       "      <td>...</td>\n",
       "      <td>Amino Acid</td>\n",
       "      <td>Phenylalanine and Tyrosine Metabolism</td>\n",
       "      <td>537-55-3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs56672945</td>\n",
       "      <td>2</td>\n",
       "      <td>ALMS1</td>\n",
       "      <td>rs56672945</td>\n",
       "      <td>Ruminococcaceae_UCG-002</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1632</th>\n",
       "      <td>metabolite_340</td>\n",
       "      <td>Ruminococcaceae|UBA1819</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>0.158704</td>\n",
       "      <td>0.024361</td>\n",
       "      <td>0.000001</td>\n",
       "      <td>0.149686</td>\n",
       "      <td>0.013943</td>\n",
       "      <td>0.163629</td>\n",
       "      <td>755</td>\n",
       "      <td>...</td>\n",
       "      <td>Amino Acid</td>\n",
       "      <td>Glycine, Serine and Threonine Metabolism</td>\n",
       "      <td>56-40-6</td>\n",
       "      <td>C00037</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>2</td>\n",
       "      <td>CPS1</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>UBA1819</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2528</th>\n",
       "      <td>metabolite_803</td>\n",
       "      <td>Ruminococcaceae|Faecalibacterium</td>\n",
       "      <td>rs1260326</td>\n",
       "      <td>-0.041809</td>\n",
       "      <td>0.02759</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>0.049705</td>\n",
       "      <td>0.014559</td>\n",
       "      <td>0.064264</td>\n",
       "      <td>247</td>\n",
       "      <td>...</td>\n",
       "      <td>Carbohydrate</td>\n",
       "      <td>Fructose, Mannose and Galactose Metabolism</td>\n",
       "      <td>3458-28-4</td>\n",
       "      <td>C00159</td>\n",
       "      <td>rs1260326</td>\n",
       "      <td>2</td>\n",
       "      <td>GCKR</td>\n",
       "      <td>rs1260326</td>\n",
       "      <td>Faecalibacterium</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18384</th>\n",
       "      <td>metabolite_100010896</td>\n",
       "      <td>Lachnospiraceae|Coprococcus_1</td>\n",
       "      <td>rs7048932</td>\n",
       "      <td>0.249172</td>\n",
       "      <td>0.021118</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>0.240782</td>\n",
       "      <td>0.011562</td>\n",
       "      <td>0.252344</td>\n",
       "      <td>607</td>\n",
       "      <td>...</td>\n",
       "      <td>Nucleotide</td>\n",
       "      <td>Pyrimidine Metabolism, Uracil containing</td>\n",
       "      <td>2140-76-3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs7048932</td>\n",
       "      <td>9</td>\n",
       "      <td>PHYHD1</td>\n",
       "      <td>rs7048932</td>\n",
       "      <td>Coprococcus_1</td>\n",
       "      <td>Lachnospiraceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1880</th>\n",
       "      <td>metabolite_100001294</td>\n",
       "      <td>Ruminococcaceae|UBA1819</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>0.190981</td>\n",
       "      <td>0.032918</td>\n",
       "      <td>0.000007</td>\n",
       "      <td>0.101733</td>\n",
       "      <td>0.012099</td>\n",
       "      <td>0.113832</td>\n",
       "      <td>755</td>\n",
       "      <td>...</td>\n",
       "      <td>Peptide</td>\n",
       "      <td>Gamma-glutamyl Amino Acid</td>\n",
       "      <td>1948-29-4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>2</td>\n",
       "      <td>CPS1</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>UBA1819</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>909</th>\n",
       "      <td>metabolite_244</td>\n",
       "      <td>Bifidobacteriaceae|Bifidobacterium</td>\n",
       "      <td>rs6800284</td>\n",
       "      <td>-0.08661</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.999882</td>\n",
       "      <td>0.047732</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.047732</td>\n",
       "      <td>477</td>\n",
       "      <td>...</td>\n",
       "      <td>Nucleotide</td>\n",
       "      <td>Pyrimidine Metabolism, Uracil containing</td>\n",
       "      <td>56-41-7;107-95-9</td>\n",
       "      <td>C00099</td>\n",
       "      <td>rs6800284</td>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs6800284</td>\n",
       "      <td>Bifidobacterium</td>\n",
       "      <td>Bifidobacteriaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22766</th>\n",
       "      <td>metabolite_999946997</td>\n",
       "      <td>Ruminococcaceae|Candidatus_Soleaferrea</td>\n",
       "      <td>rs765285</td>\n",
       "      <td>0.161182</td>\n",
       "      <td>-0.000002</td>\n",
       "      <td>0.999889</td>\n",
       "      <td>0.052567</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.052567</td>\n",
       "      <td>324</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs765285</td>\n",
       "      <td>6</td>\n",
       "      <td>SLC17A1</td>\n",
       "      <td>rs765285</td>\n",
       "      <td>Candidatus_Soleaferrea</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14783</th>\n",
       "      <td>metabolite_100002259</td>\n",
       "      <td>Lachnospiraceae|Fusicatenibacter</td>\n",
       "      <td>rs12134854</td>\n",
       "      <td>-0.182032</td>\n",
       "      <td>-0.000001</td>\n",
       "      <td>0.999923</td>\n",
       "      <td>0.062553</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.062553</td>\n",
       "      <td>819</td>\n",
       "      <td>...</td>\n",
       "      <td>Lipid</td>\n",
       "      <td>Fatty Acid Metabolism(Acyl Carnitine)</td>\n",
       "      <td>98930-66-6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs12134854</td>\n",
       "      <td>1</td>\n",
       "      <td>SLC44A5</td>\n",
       "      <td>rs12134854</td>\n",
       "      <td>Fusicatenibacter</td>\n",
       "      <td>Lachnospiraceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16656</th>\n",
       "      <td>metabolite_100009078</td>\n",
       "      <td>Erysipelotrichaceae|Erysipelotrichaceae_UCG-003</td>\n",
       "      <td>rs35853021</td>\n",
       "      <td>0.156265</td>\n",
       "      <td>-0.000001</td>\n",
       "      <td>0.999923</td>\n",
       "      <td>0.04425</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.04425</td>\n",
       "      <td>638</td>\n",
       "      <td>...</td>\n",
       "      <td>Lipid</td>\n",
       "      <td>Phospholipid Metabolism</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs35853021</td>\n",
       "      <td>15</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs35853021</td>\n",
       "      <td>Erysipelotrichaceae_UCG-003</td>\n",
       "      <td>Erysipelotrichaceae</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4771</th>\n",
       "      <td>metabolite_2054</td>\n",
       "      <td>Ruminococcaceae|Subdoligranulum</td>\n",
       "      <td>rs1799958</td>\n",
       "      <td>0.538413</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.485331</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>0.485331</td>\n",
       "      <td>857</td>\n",
       "      <td>...</td>\n",
       "      <td>Amino Acid</td>\n",
       "      <td>Leucine, Isoleucine and Valine Metabolism</td>\n",
       "      <td>601-75-2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>rs1799958</td>\n",
       "      <td>12</td>\n",
       "      <td>ACADS</td>\n",
       "      <td>rs1799958</td>\n",
       "      <td>Subdoligranulum</td>\n",
       "      <td>Ruminococcaceae</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>23126 rows × 25 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 metabolite                                            taxon  \\\n",
       "9687   metabolite_100001104          Ruminococcaceae|Ruminococcaceae_UCG-002   \n",
       "1632         metabolite_340                          Ruminococcaceae|UBA1819   \n",
       "2528         metabolite_803                 Ruminococcaceae|Faecalibacterium   \n",
       "18384  metabolite_100010896                    Lachnospiraceae|Coprococcus_1   \n",
       "1880   metabolite_100001294                          Ruminococcaceae|UBA1819   \n",
       "...                     ...                                              ...   \n",
       "909          metabolite_244               Bifidobacteriaceae|Bifidobacterium   \n",
       "22766  metabolite_999946997           Ruminococcaceae|Candidatus_Soleaferrea   \n",
       "14783  metabolite_100002259                 Lachnospiraceae|Fusicatenibacter   \n",
       "16656  metabolite_100009078  Erysipelotrichaceae|Erysipelotrichaceae_UCG-003   \n",
       "4771        metabolite_2054                  Ruminococcaceae|Subdoligranulum   \n",
       "\n",
       "             snps  snp_beta interaction_beta         p baseline_r2  \\\n",
       "9687   rs56672945  0.309384        -0.032553  0.000001    0.189084   \n",
       "1632    rs1047891  0.158704         0.024361  0.000001    0.149686   \n",
       "2528    rs1260326 -0.041809          0.02759  0.000002    0.049705   \n",
       "18384   rs7048932  0.249172         0.021118  0.000004    0.240782   \n",
       "1880    rs1047891  0.190981         0.032918  0.000007    0.101733   \n",
       "...           ...       ...              ...       ...         ...   \n",
       "909     rs6800284  -0.08661             -0.0  0.999882    0.047732   \n",
       "22766    rs765285  0.161182        -0.000002  0.999889    0.052567   \n",
       "14783  rs12134854 -0.182032        -0.000001  0.999923    0.062553   \n",
       "16656  rs35853021  0.156265        -0.000001  0.999923     0.04425   \n",
       "4771    rs1799958  0.538413             -0.0       1.0    0.485331   \n",
       "\n",
       "      interaction_r2   full_r2 n_major  ... SUPER_PATHWAY  \\\n",
       "9687        0.014297  0.203381     944  ...    Amino Acid   \n",
       "1632        0.013943  0.163629     755  ...    Amino Acid   \n",
       "2528        0.014559  0.064264     247  ...  Carbohydrate   \n",
       "18384       0.011562  0.252344     607  ...    Nucleotide   \n",
       "1880        0.012099  0.113832     755  ...       Peptide   \n",
       "...              ...       ...     ...  ...           ...   \n",
       "909              0.0  0.047732     477  ...    Nucleotide   \n",
       "22766            0.0  0.052567     324  ...           NaN   \n",
       "14783            0.0  0.062553     819  ...         Lipid   \n",
       "16656            0.0   0.04425     638  ...         Lipid   \n",
       "4771            -0.0  0.485331     857  ...    Amino Acid   \n",
       "\n",
       "                                      SUB_PATHWAY               CAS    KEGG  \\\n",
       "9687        Phenylalanine and Tyrosine Metabolism          537-55-3     NaN   \n",
       "1632     Glycine, Serine and Threonine Metabolism           56-40-6  C00037   \n",
       "2528   Fructose, Mannose and Galactose Metabolism         3458-28-4  C00159   \n",
       "18384    Pyrimidine Metabolism, Uracil containing         2140-76-3     NaN   \n",
       "1880                    Gamma-glutamyl Amino Acid         1948-29-4     NaN   \n",
       "...                                           ...               ...     ...   \n",
       "909      Pyrimidine Metabolism, Uracil containing  56-41-7;107-95-9  C00099   \n",
       "22766                                         NaN               NaN     NaN   \n",
       "14783       Fatty Acid Metabolism(Acyl Carnitine)        98930-66-6     NaN   \n",
       "16656                     Phospholipid Metabolism               NaN     NaN   \n",
       "4771    Leucine, Isoleucine and Valine Metabolism          601-75-2     NaN   \n",
       "\n",
       "              SNP CHR    genes        rsid                        genus  \\\n",
       "9687   rs56672945   2    ALMS1  rs56672945      Ruminococcaceae_UCG-002   \n",
       "1632    rs1047891   2     CPS1   rs1047891                      UBA1819   \n",
       "2528    rs1260326   2     GCKR   rs1260326             Faecalibacterium   \n",
       "18384   rs7048932   9   PHYHD1   rs7048932                Coprococcus_1   \n",
       "1880    rs1047891   2     CPS1   rs1047891                      UBA1819   \n",
       "...           ...  ..      ...         ...                          ...   \n",
       "909     rs6800284   3      NaN   rs6800284              Bifidobacterium   \n",
       "22766    rs765285   6  SLC17A1    rs765285       Candidatus_Soleaferrea   \n",
       "14783  rs12134854   1  SLC44A5  rs12134854             Fusicatenibacter   \n",
       "16656  rs35853021  15      NaN  rs35853021  Erysipelotrichaceae_UCG-003   \n",
       "4771    rs1799958  12    ACADS   rs1799958              Subdoligranulum   \n",
       "\n",
       "                    family  \n",
       "9687       Ruminococcaceae  \n",
       "1632       Ruminococcaceae  \n",
       "2528       Ruminococcaceae  \n",
       "18384      Lachnospiraceae  \n",
       "1880       Ruminococcaceae  \n",
       "...                    ...  \n",
       "909     Bifidobacteriaceae  \n",
       "22766      Ruminococcaceae  \n",
       "14783      Lachnospiraceae  \n",
       "16656  Erysipelotrichaceae  \n",
       "4771       Ruminococcaceae  \n",
       "\n",
       "[23126 rows x 25 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metabolite_meta = adi.get_snapshot(\"metabolomics_metadata\", clean=True) [[\"CHEMICAL_ID\", \"BIOCHEMICAL_NAME\", \"SUPER_PATHWAY\", \"SUB_PATHWAY\", \"CAS\", \"KEGG\"]]\n",
    "metabolite_meta[\"metabolite\"] = \"metabolite_\" + metabolite_meta.CHEMICAL_ID.astype(str)\n",
    "merged = pd.merge(results, metabolite_meta, on=\"metabolite\")\n",
    "merged = pd.merge(merged, sig_snps[[\"SNP\", \"CHR\", \"genes\", \"rsid\"]].drop_duplicates(), left_on=\"snps\", right_on=\"SNP\")\n",
    "merged[\"genus\"] = merged.taxon.str.split(\"|\").str[1]\n",
    "merged[\"family\"] = merged.taxon.str.split(\"|\").str[0]\n",
    "\n",
    "merged.sort_values(by=[\"p\", \"interaction_r2\"], inplace=True)\n",
    "merged.to_csv(\"data/interaction_results_annotated.csv\", index=False)\n",
    "merged"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hybrid-fitting",
   "metadata": {},
   "source": [
    "Let's have a look at the significant interactions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "coordinate-market",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected = merged[merged.q < 0.05].sort_values(by=\"interaction_r2\", ascending=False)[[\n",
    "    \"BIOCHEMICAL_NAME\", \"interaction_r2\", \"interaction_beta\", \"baseline_r2\", \n",
    "    \"SNP\", \"genes\", \"n_major\", \"n_heterozygous\", \"n_homozygous\", \"taxon\", \"p\", \"q\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "iraqi-russia",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ruminococcaceae|Ruminococcaceae_UCG-013     20\n",
       "Ruminococcaceae|Angelakisella               13\n",
       "Ruminococcaceae|UBA1819                     13\n",
       "Ruminococcaceae|DTU089                      12\n",
       "Veillonellaceae|Veillonella                 12\n",
       "                                            ..\n",
       "Ruminococcaceae|Caproiciproducens            1\n",
       "Desulfovibrionaceae|Bilophila                1\n",
       "Eggerthellaceae|Adlercreutzia                1\n",
       "Acidaminococcaceae|Phascolarctobacterium     1\n",
       "Ruminococcaceae|Ruminococcus_2               1\n",
       "Name: taxon, Length: 79, dtype: int64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected.taxon.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "turkish-execution",
   "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>BIOCHEMICAL_NAME</th>\n",
       "      <th>interaction_r2</th>\n",
       "      <th>interaction_beta</th>\n",
       "      <th>baseline_r2</th>\n",
       "      <th>SNP</th>\n",
       "      <th>genes</th>\n",
       "      <th>n_major</th>\n",
       "      <th>n_heterozygous</th>\n",
       "      <th>n_homozygous</th>\n",
       "      <th>taxon</th>\n",
       "      <th>p</th>\n",
       "      <th>q</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2528</th>\n",
       "      <td>mannose</td>\n",
       "      <td>0.014559</td>\n",
       "      <td>0.02759</td>\n",
       "      <td>0.049705</td>\n",
       "      <td>rs1260326</td>\n",
       "      <td>GCKR</td>\n",
       "      <td>247</td>\n",
       "      <td>758</td>\n",
       "      <td>564</td>\n",
       "      <td>Ruminococcaceae|Faecalibacterium</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>0.000003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9687</th>\n",
       "      <td>N-acetyltyrosine</td>\n",
       "      <td>0.014297</td>\n",
       "      <td>-0.032553</td>\n",
       "      <td>0.189084</td>\n",
       "      <td>rs56672945</td>\n",
       "      <td>ALMS1</td>\n",
       "      <td>944</td>\n",
       "      <td>529</td>\n",
       "      <td>96</td>\n",
       "      <td>Ruminococcaceae|Ruminococcaceae_UCG-002</td>\n",
       "      <td>0.000001</td>\n",
       "      <td>0.000001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1632</th>\n",
       "      <td>glycine</td>\n",
       "      <td>0.013943</td>\n",
       "      <td>0.024361</td>\n",
       "      <td>0.149686</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>CPS1</td>\n",
       "      <td>755</td>\n",
       "      <td>663</td>\n",
       "      <td>151</td>\n",
       "      <td>Ruminococcaceae|UBA1819</td>\n",
       "      <td>0.000001</td>\n",
       "      <td>0.000001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3601</th>\n",
       "      <td>cys-gly, oxidized</td>\n",
       "      <td>0.012553</td>\n",
       "      <td>0.041491</td>\n",
       "      <td>0.052908</td>\n",
       "      <td>rs258341</td>\n",
       "      <td>DPEP1</td>\n",
       "      <td>261</td>\n",
       "      <td>785</td>\n",
       "      <td>523</td>\n",
       "      <td>Ruminococcaceae|Faecalibacterium</td>\n",
       "      <td>0.000008</td>\n",
       "      <td>0.000016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1880</th>\n",
       "      <td>gamma-glutamylglycine</td>\n",
       "      <td>0.012099</td>\n",
       "      <td>0.032918</td>\n",
       "      <td>0.101733</td>\n",
       "      <td>rs1047891</td>\n",
       "      <td>CPS1</td>\n",
       "      <td>755</td>\n",
       "      <td>663</td>\n",
       "      <td>151</td>\n",
       "      <td>Ruminococcaceae|UBA1819</td>\n",
       "      <td>0.000007</td>\n",
       "      <td>0.000009</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           BIOCHEMICAL_NAME interaction_r2 interaction_beta baseline_r2  \\\n",
       "2528                mannose       0.014559          0.02759    0.049705   \n",
       "9687       N-acetyltyrosine       0.014297        -0.032553    0.189084   \n",
       "1632                glycine       0.013943         0.024361    0.149686   \n",
       "3601      cys-gly, oxidized       0.012553         0.041491    0.052908   \n",
       "1880  gamma-glutamylglycine       0.012099         0.032918    0.101733   \n",
       "\n",
       "             SNP  genes n_major n_heterozygous n_homozygous  \\\n",
       "2528   rs1260326   GCKR     247            758          564   \n",
       "9687  rs56672945  ALMS1     944            529           96   \n",
       "1632   rs1047891   CPS1     755            663          151   \n",
       "3601    rs258341  DPEP1     261            785          523   \n",
       "1880   rs1047891   CPS1     755            663          151   \n",
       "\n",
       "                                        taxon         p         q  \n",
       "2528         Ruminococcaceae|Faecalibacterium  0.000002  0.000003  \n",
       "9687  Ruminococcaceae|Ruminococcaceae_UCG-002  0.000001  0.000001  \n",
       "1632                  Ruminococcaceae|UBA1819  0.000001  0.000001  \n",
       "3601         Ruminococcaceae|Faecalibacterium  0.000008  0.000016  \n",
       "1880                  Ruminococcaceae|UBA1819  0.000007  0.000009  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "desirable-emergency",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(466, 12)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "suited-exploration",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected.BIOCHEMICAL_NAME.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "environmental-frontier",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.004516474711522251"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected.interaction_r2.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "changed-import",
   "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}