microbe_associations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "sized-lawyer",
   "metadata": {},
   "source": [
    "# Identify bacterial genera associated with metabolites\n",
    "\n",
    "We will now select the microbes significantly associated with each metabolite. We will use the genus level here since that is usually the lowest taxonomic rank we can confidently identify with 16S amplicon sequencing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "loved-boulder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "410"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "train = pd.read_csv(\"data/train.csv\")\n",
    "valid = pd.read_csv(\"data/valid.csv\")\n",
    "abundance = pd.read_csv(\"/proj/arivale/microbiome/16S_processed/asvs.csv\")\n",
    "abundance[\"stool_sample_id\"] = abundance.id.str.split(\"|\").str[0]\n",
    "abundance = abundance[abundance.stool_sample_id.isin(train.stool_sample_id) | abundance.stool_sample_id.isin(valid.stool_sample_id)]\n",
    "\n",
    "taxonomy = pd.read_csv(\"/proj/arivale/microbiome/16S_processed/taxonomy.csv\")\n",
    "taxa_compact = taxonomy.Family + \"|\" + taxonomy.Genus\n",
    "taxa_compact.index = taxonomy.id\n",
    "abundance[\"taxon\"] = taxa_compact[abundance.hash].values\n",
    "genera = abundance.groupby([\"stool_sample_id\", \"taxon\"]).sum().reset_index()\n",
    "genera.taxon.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "nominated-awareness",
   "metadata": {},
   "source": [
    "We also filter out genera that were not observed in at least 50% of all individuals and that have less than 10 detected reads per sample on average. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "precise-sharing",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "82"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prevalence = genera.groupby(\"taxon\").apply(lambda df: sum(df[\"count\"] > 0)) / genera.stool_sample_id.nunique()\n",
    "means = genera.groupby(\"taxon\")[\"count\"].mean()\n",
    "good = (prevalence > 0.5) & (means > 10)\n",
    "genera_filtered = genera[genera.taxon.isin(good.index[good])]\n",
    "genera_filtered.taxon.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "forward-crisis",
   "metadata": {},
   "source": [
    "Now we will CLR transform the abundances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "chubby-advertising",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1569, 82)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def clr(x, pseudo):\n",
    "    if pseudo == 0:\n",
    "        out = x.copy()\n",
    "        out[x == 0] = float(\"nan\")\n",
    "        out[x > 0] = np.log(x[x > 0]) - np.mean(np.log(x[x > 0])) \n",
    "    else:\n",
    "        out = np.log(x + pseudo) - np.mean(np.log(x + pseudo))\n",
    "    return out\n",
    "\n",
    "genera_wide = pd.pivot_table(genera_filtered, index=\"stool_sample_id\", columns=\"taxon\", values=\"count\", fill_value=0)\n",
    "genera_wide_clr = genera_wide.apply(lambda x: clr(x, 1), axis=1)\n",
    "genera_wide_clr.to_csv(\"data/genera_clr_filtered.csv\")\n",
    "\n",
    "genera_wide_clr_drop = genera_wide.apply(lambda x: clr(x, 0), axis=1)\n",
    "genera_wide_clr_drop.to_csv(\"data/genera_clr_filtered_nonzero.csv\")\n",
    "genera_wide_clr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "chinese-lobby",
   "metadata": {},
   "outputs": [],
   "source": [
    "metabolites_and_microbes = pd.merge(train, genera_wide_clr, left_on=\"stool_sample_id\", right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "swedish-mumbai",
   "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>public_client_id</th>\n",
       "      <th>genome_id</th>\n",
       "      <th>blood_sample_id</th>\n",
       "      <th>blood_days_in_program</th>\n",
       "      <th>stool_sample_id</th>\n",
       "      <th>stool_days_in_program</th>\n",
       "      <th>sex</th>\n",
       "      <th>age</th>\n",
       "      <th>stool_vendor</th>\n",
       "      <th>WEIGHT_CALC</th>\n",
       "      <th>...</th>\n",
       "      <th>Ruminococcaceae|Ruminococcaceae_UCG-014</th>\n",
       "      <th>Ruminococcaceae|Ruminococcus_1</th>\n",
       "      <th>Ruminococcaceae|Ruminococcus_2</th>\n",
       "      <th>Ruminococcaceae|Subdoligranulum</th>\n",
       "      <th>Ruminococcaceae|UBA1819</th>\n",
       "      <th>Streptococcaceae|Lactococcus</th>\n",
       "      <th>Streptococcaceae|Streptococcus</th>\n",
       "      <th>Tannerellaceae|Parabacteroides</th>\n",
       "      <th>Veillonellaceae|Dialister</th>\n",
       "      <th>Veillonellaceae|Veillonella</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1241633</td>\n",
       "      <td>A789AU958-007</td>\n",
       "      <td>A789AU992-002</td>\n",
       "      <td>36.0</td>\n",
       "      <td>AV15-1655</td>\n",
       "      <td>36.0</td>\n",
       "      <td>M</td>\n",
       "      <td>71.0</td>\n",
       "      <td>SecondGenome</td>\n",
       "      <td>179.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.462865</td>\n",
       "      <td>3.690187</td>\n",
       "      <td>-0.061668</td>\n",
       "      <td>5.699964</td>\n",
       "      <td>1.513869</td>\n",
       "      <td>-1.265641</td>\n",
       "      <td>2.477306</td>\n",
       "      <td>1.506948</td>\n",
       "      <td>4.792963</td>\n",
       "      <td>-0.284811</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1020175</td>\n",
       "      <td>DS-268800</td>\n",
       "      <td>A430BH742-003</td>\n",
       "      <td>8.0</td>\n",
       "      <td>22001612561075</td>\n",
       "      <td>13.0</td>\n",
       "      <td>F</td>\n",
       "      <td>41.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>192.8</td>\n",
       "      <td>...</td>\n",
       "      <td>-4.151371</td>\n",
       "      <td>2.642095</td>\n",
       "      <td>-4.151371</td>\n",
       "      <td>2.343894</td>\n",
       "      <td>0.892054</td>\n",
       "      <td>-0.413702</td>\n",
       "      <td>0.825362</td>\n",
       "      <td>5.826435</td>\n",
       "      <td>3.427797</td>\n",
       "      <td>-1.015877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1230573</td>\n",
       "      <td>DS-284622</td>\n",
       "      <td>A773BJ600-002</td>\n",
       "      <td>4.0</td>\n",
       "      <td>22001701512104</td>\n",
       "      <td>22.0</td>\n",
       "      <td>F</td>\n",
       "      <td>60.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>152.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.840533</td>\n",
       "      <td>-4.836265</td>\n",
       "      <td>-4.836265</td>\n",
       "      <td>-4.836265</td>\n",
       "      <td>0.436734</td>\n",
       "      <td>-1.122693</td>\n",
       "      <td>-1.029603</td>\n",
       "      <td>3.860078</td>\n",
       "      <td>3.405438</td>\n",
       "      <td>-4.836265</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1687802</td>\n",
       "      <td>A719AT743-004</td>\n",
       "      <td>A719AT785-002</td>\n",
       "      <td>88.0</td>\n",
       "      <td>AV15-1320</td>\n",
       "      <td>88.0</td>\n",
       "      <td>F</td>\n",
       "      <td>53.0</td>\n",
       "      <td>SecondGenome</td>\n",
       "      <td>242.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.409260</td>\n",
       "      <td>3.829955</td>\n",
       "      <td>3.229308</td>\n",
       "      <td>1.300270</td>\n",
       "      <td>0.598073</td>\n",
       "      <td>1.155088</td>\n",
       "      <td>2.365292</td>\n",
       "      <td>-3.409260</td>\n",
       "      <td>-3.409260</td>\n",
       "      <td>0.117101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1864964</td>\n",
       "      <td>DS-267221</td>\n",
       "      <td>A862BF047-008</td>\n",
       "      <td>14.0</td>\n",
       "      <td>22001612561843</td>\n",
       "      <td>15.0</td>\n",
       "      <td>F</td>\n",
       "      <td>43.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.535311</td>\n",
       "      <td>4.653378</td>\n",
       "      <td>1.226863</td>\n",
       "      <td>2.005952</td>\n",
       "      <td>-0.762722</td>\n",
       "      <td>-3.535311</td>\n",
       "      <td>0.942026</td>\n",
       "      <td>4.230682</td>\n",
       "      <td>-3.535311</td>\n",
       "      <td>-0.827261</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>995</th>\n",
       "      <td>1035670</td>\n",
       "      <td>DS-269430</td>\n",
       "      <td>A808BG205-003</td>\n",
       "      <td>33.0</td>\n",
       "      <td>22001612562579</td>\n",
       "      <td>9.0</td>\n",
       "      <td>F</td>\n",
       "      <td>33.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>155.2</td>\n",
       "      <td>...</td>\n",
       "      <td>-5.483529</td>\n",
       "      <td>4.261138</td>\n",
       "      <td>3.655530</td>\n",
       "      <td>4.293317</td>\n",
       "      <td>-0.608332</td>\n",
       "      <td>-2.392487</td>\n",
       "      <td>0.461891</td>\n",
       "      <td>1.905417</td>\n",
       "      <td>2.450267</td>\n",
       "      <td>-2.187693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>1013950</td>\n",
       "      <td>A752AT594-004</td>\n",
       "      <td>A752AT607-002</td>\n",
       "      <td>77.0</td>\n",
       "      <td>AV15-1319</td>\n",
       "      <td>77.0</td>\n",
       "      <td>F</td>\n",
       "      <td>42.0</td>\n",
       "      <td>SecondGenome</td>\n",
       "      <td>126.0</td>\n",
       "      <td>...</td>\n",
       "      <td>3.016214</td>\n",
       "      <td>2.118971</td>\n",
       "      <td>4.178741</td>\n",
       "      <td>3.622815</td>\n",
       "      <td>0.170962</td>\n",
       "      <td>-5.147158</td>\n",
       "      <td>2.745668</td>\n",
       "      <td>1.831056</td>\n",
       "      <td>-5.147158</td>\n",
       "      <td>-5.147158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>1244110</td>\n",
       "      <td>A939AV031-007</td>\n",
       "      <td>A473AS982-002</td>\n",
       "      <td>52.0</td>\n",
       "      <td>AV15-1143</td>\n",
       "      <td>52.0</td>\n",
       "      <td>F</td>\n",
       "      <td>62.0</td>\n",
       "      <td>SecondGenome</td>\n",
       "      <td>126.0</td>\n",
       "      <td>...</td>\n",
       "      <td>2.807030</td>\n",
       "      <td>1.858001</td>\n",
       "      <td>2.510743</td>\n",
       "      <td>1.855598</td>\n",
       "      <td>-2.277166</td>\n",
       "      <td>-1.717550</td>\n",
       "      <td>1.581457</td>\n",
       "      <td>1.387677</td>\n",
       "      <td>2.798008</td>\n",
       "      <td>-5.272898</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>1329753</td>\n",
       "      <td>DS-265707</td>\n",
       "      <td>A579BF997-003</td>\n",
       "      <td>3.0</td>\n",
       "      <td>22001612560465</td>\n",
       "      <td>15.0</td>\n",
       "      <td>M</td>\n",
       "      <td>40.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>269.0</td>\n",
       "      <td>...</td>\n",
       "      <td>2.352577</td>\n",
       "      <td>3.215964</td>\n",
       "      <td>3.505026</td>\n",
       "      <td>3.807427</td>\n",
       "      <td>-3.081145</td>\n",
       "      <td>-3.081145</td>\n",
       "      <td>0.870099</td>\n",
       "      <td>1.429714</td>\n",
       "      <td>-3.081145</td>\n",
       "      <td>-3.081145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>999</th>\n",
       "      <td>1028150</td>\n",
       "      <td>DS-268107</td>\n",
       "      <td>A202BI892-020</td>\n",
       "      <td>30.0</td>\n",
       "      <td>22001612562223</td>\n",
       "      <td>25.0</td>\n",
       "      <td>M</td>\n",
       "      <td>57.0</td>\n",
       "      <td>DNAGenotek</td>\n",
       "      <td>225.0</td>\n",
       "      <td>...</td>\n",
       "      <td>2.794583</td>\n",
       "      <td>4.623733</td>\n",
       "      <td>3.963858</td>\n",
       "      <td>-0.579437</td>\n",
       "      <td>0.609202</td>\n",
       "      <td>-0.051369</td>\n",
       "      <td>1.912690</td>\n",
       "      <td>3.838356</td>\n",
       "      <td>3.229638</td>\n",
       "      <td>-4.408078</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1000 rows × 1030 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     public_client_id      genome_id blood_sample_id  blood_days_in_program  \\\n",
       "0             1241633  A789AU958-007   A789AU992-002                   36.0   \n",
       "1             1020175      DS-268800   A430BH742-003                    8.0   \n",
       "2             1230573      DS-284622   A773BJ600-002                    4.0   \n",
       "3             1687802  A719AT743-004   A719AT785-002                   88.0   \n",
       "4             1864964      DS-267221   A862BF047-008                   14.0   \n",
       "..                ...            ...             ...                    ...   \n",
       "995           1035670      DS-269430   A808BG205-003                   33.0   \n",
       "996           1013950  A752AT594-004   A752AT607-002                   77.0   \n",
       "997           1244110  A939AV031-007   A473AS982-002                   52.0   \n",
       "998           1329753      DS-265707   A579BF997-003                    3.0   \n",
       "999           1028150      DS-268107   A202BI892-020                   30.0   \n",
       "\n",
       "    stool_sample_id  stool_days_in_program sex   age  stool_vendor  \\\n",
       "0         AV15-1655                   36.0   M  71.0  SecondGenome   \n",
       "1    22001612561075                   13.0   F  41.0    DNAGenotek   \n",
       "2    22001701512104                   22.0   F  60.0    DNAGenotek   \n",
       "3         AV15-1320                   88.0   F  53.0  SecondGenome   \n",
       "4    22001612561843                   15.0   F  43.0    DNAGenotek   \n",
       "..              ...                    ...  ..   ...           ...   \n",
       "995  22001612562579                    9.0   F  33.0    DNAGenotek   \n",
       "996       AV15-1319                   77.0   F  42.0  SecondGenome   \n",
       "997       AV15-1143                   52.0   F  62.0  SecondGenome   \n",
       "998  22001612560465                   15.0   M  40.0    DNAGenotek   \n",
       "999  22001612562223                   25.0   M  57.0    DNAGenotek   \n",
       "\n",
       "     WEIGHT_CALC  ...  Ruminococcaceae|Ruminococcaceae_UCG-014  \\\n",
       "0          179.0  ...                                -3.462865   \n",
       "1          192.8  ...                                -4.151371   \n",
       "2          152.0  ...                                -1.840533   \n",
       "3          242.0  ...                                -3.409260   \n",
       "4          170.0  ...                                -3.535311   \n",
       "..           ...  ...                                      ...   \n",
       "995        155.2  ...                                -5.483529   \n",
       "996        126.0  ...                                 3.016214   \n",
       "997        126.0  ...                                 2.807030   \n",
       "998        269.0  ...                                 2.352577   \n",
       "999        225.0  ...                                 2.794583   \n",
       "\n",
       "    Ruminococcaceae|Ruminococcus_1 Ruminococcaceae|Ruminococcus_2  \\\n",
       "0                         3.690187                      -0.061668   \n",
       "1                         2.642095                      -4.151371   \n",
       "2                        -4.836265                      -4.836265   \n",
       "3                         3.829955                       3.229308   \n",
       "4                         4.653378                       1.226863   \n",
       "..                             ...                            ...   \n",
       "995                       4.261138                       3.655530   \n",
       "996                       2.118971                       4.178741   \n",
       "997                       1.858001                       2.510743   \n",
       "998                       3.215964                       3.505026   \n",
       "999                       4.623733                       3.963858   \n",
       "\n",
       "     Ruminococcaceae|Subdoligranulum  Ruminococcaceae|UBA1819  \\\n",
       "0                           5.699964                 1.513869   \n",
       "1                           2.343894                 0.892054   \n",
       "2                          -4.836265                 0.436734   \n",
       "3                           1.300270                 0.598073   \n",
       "4                           2.005952                -0.762722   \n",
       "..                               ...                      ...   \n",
       "995                         4.293317                -0.608332   \n",
       "996                         3.622815                 0.170962   \n",
       "997                         1.855598                -2.277166   \n",
       "998                         3.807427                -3.081145   \n",
       "999                        -0.579437                 0.609202   \n",
       "\n",
       "     Streptococcaceae|Lactococcus  Streptococcaceae|Streptococcus  \\\n",
       "0                       -1.265641                        2.477306   \n",
       "1                       -0.413702                        0.825362   \n",
       "2                       -1.122693                       -1.029603   \n",
       "3                        1.155088                        2.365292   \n",
       "4                       -3.535311                        0.942026   \n",
       "..                            ...                             ...   \n",
       "995                     -2.392487                        0.461891   \n",
       "996                     -5.147158                        2.745668   \n",
       "997                     -1.717550                        1.581457   \n",
       "998                     -3.081145                        0.870099   \n",
       "999                     -0.051369                        1.912690   \n",
       "\n",
       "     Tannerellaceae|Parabacteroides  Veillonellaceae|Dialister  \\\n",
       "0                          1.506948                   4.792963   \n",
       "1                          5.826435                   3.427797   \n",
       "2                          3.860078                   3.405438   \n",
       "3                         -3.409260                  -3.409260   \n",
       "4                          4.230682                  -3.535311   \n",
       "..                              ...                        ...   \n",
       "995                        1.905417                   2.450267   \n",
       "996                        1.831056                  -5.147158   \n",
       "997                        1.387677                   2.798008   \n",
       "998                        1.429714                  -3.081145   \n",
       "999                        3.838356                   3.229638   \n",
       "\n",
       "     Veillonellaceae|Veillonella  \n",
       "0                      -0.284811  \n",
       "1                      -1.015877  \n",
       "2                      -4.836265  \n",
       "3                       0.117101  \n",
       "4                      -0.827261  \n",
       "..                           ...  \n",
       "995                    -2.187693  \n",
       "996                    -5.147158  \n",
       "997                    -5.272898  \n",
       "998                    -3.081145  \n",
       "999                    -4.408078  \n",
       "\n",
       "[1000 rows x 1030 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metabolites_and_microbes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "civil-worth",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37e4b3d1d01241c68880b48a77e4b64c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from itertools import product\n",
    "from rich.progress import track\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "mets = metabolites_and_microbes.columns[metabolites_and_microbes.columns.str.startswith(\"metabolite\")]\n",
    "mics = genera_wide_clr.columns\n",
    "\n",
    "combinations = list(product(mets, mics))\n",
    "\n",
    "metrics = pd.DataFrame(columns=[\"r\", \"p\", \"n\"], dtype=\"float64\", index=combinations)\n",
    "formula = \"{} ~ Q('{}')\"\n",
    "for comb in track(combinations):\n",
    "    data = metabolites_and_microbes[list(comb)].dropna()\n",
    "    ptest = pearsonr(data.iloc[:, 0], data.iloc[:, 1])\n",
    "    metrics.loc[[comb], \"r\"] = ptest[0]\n",
    "    metrics.loc[[comb], \"p\"] = ptest[1]\n",
    "    metrics.loc[[comb], \"n\"] = data.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "engaging-stability",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.multitest import fdrcorrection\n",
    "\n",
    "metrics[\"q\"] = fdrcorrection(metrics.p)[1]\n",
    "metrics[\"r2\"] = metrics.r.pow(2)\n",
    "metrics[\"metabolite\"] = metrics.index.map(lambda x: x[0])\n",
    "metrics[\"taxon\"] = metrics.index.map(lambda x: x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "exact-mexico",
   "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>r</th>\n",
       "      <th>p</th>\n",
       "      <th>n</th>\n",
       "      <th>q</th>\n",
       "      <th>r2</th>\n",
       "      <th>metabolite</th>\n",
       "      <th>taxon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(metabolite_100002253, Family_XIII|Family_XIII_UCG-001)</th>\n",
       "      <td>0.508087</td>\n",
       "      <td>5.617983e-60</td>\n",
       "      <td>896.0</td>\n",
       "      <td>4.284274e-55</td>\n",
       "      <td>0.258152</td>\n",
       "      <td>metabolite_100002253</td>\n",
       "      <td>Family_XIII|Family_XIII_UCG-001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100002488, Ruminococcaceae|Ruminococcaceae_NK4A214_group)</th>\n",
       "      <td>-0.476256</td>\n",
       "      <td>7.957527e-54</td>\n",
       "      <td>930.0</td>\n",
       "      <td>2.817885e-49</td>\n",
       "      <td>0.226820</td>\n",
       "      <td>metabolite_100002488</td>\n",
       "      <td>Ruminococcaceae|Ruminococcaceae_NK4A214_group</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100002488, Ruminococcaceae|Flavonifractor)</th>\n",
       "      <td>0.475678</td>\n",
       "      <td>1.108531e-53</td>\n",
       "      <td>930.0</td>\n",
       "      <td>2.817885e-49</td>\n",
       "      <td>0.226270</td>\n",
       "      <td>metabolite_100002488</td>\n",
       "      <td>Ruminococcaceae|Flavonifractor</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100002488, Lachnospiraceae|Lachnoclostridium)</th>\n",
       "      <td>0.468988</td>\n",
       "      <td>4.900222e-52</td>\n",
       "      <td>930.0</td>\n",
       "      <td>9.342274e-48</td>\n",
       "      <td>0.219949</td>\n",
       "      <td>metabolite_100002488</td>\n",
       "      <td>Lachnospiraceae|Lachnoclostridium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100000010, Family_XIII|Family_XIII_UCG-001)</th>\n",
       "      <td>0.442081</td>\n",
       "      <td>4.059205e-43</td>\n",
       "      <td>874.0</td>\n",
       "      <td>6.191099e-39</td>\n",
       "      <td>0.195436</td>\n",
       "      <td>metabolite_100000010</td>\n",
       "      <td>Family_XIII|Family_XIII_UCG-001</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100020536, Lachnospiraceae|Lachnoclostridium)</th>\n",
       "      <td>-0.272563</td>\n",
       "      <td>9.977384e-17</td>\n",
       "      <td>896.0</td>\n",
       "      <td>7.925785e-14</td>\n",
       "      <td>0.074291</td>\n",
       "      <td>metabolite_100020536</td>\n",
       "      <td>Lachnospiraceae|Lachnoclostridium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100001417, Lachnospiraceae|Lachnospiraceae_UCG-004)</th>\n",
       "      <td>-0.265274</td>\n",
       "      <td>1.484930e-16</td>\n",
       "      <td>937.0</td>\n",
       "      <td>1.167431e-13</td>\n",
       "      <td>0.070370</td>\n",
       "      <td>metabolite_100001417</td>\n",
       "      <td>Lachnospiraceae|Lachnospiraceae_UCG-004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100001083, Erysipelotrichaceae|Erysipelatoclostridium)</th>\n",
       "      <td>-0.263484</td>\n",
       "      <td>1.863074e-16</td>\n",
       "      <td>944.0</td>\n",
       "      <td>1.449776e-13</td>\n",
       "      <td>0.069424</td>\n",
       "      <td>metabolite_100001083</td>\n",
       "      <td>Erysipelotrichaceae|Erysipelatoclostridium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100000014, Ruminococcaceae|Ruminococcaceae_UCG-002)</th>\n",
       "      <td>0.260556</td>\n",
       "      <td>2.574030e-16</td>\n",
       "      <td>957.0</td>\n",
       "      <td>1.982783e-13</td>\n",
       "      <td>0.067889</td>\n",
       "      <td>metabolite_100000014</td>\n",
       "      <td>Ruminococcaceae|Ruminococcaceae_UCG-002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(metabolite_100000010, Lachnospiraceae|Sellimonas)</th>\n",
       "      <td>-0.271929</td>\n",
       "      <td>2.783801e-16</td>\n",
       "      <td>874.0</td>\n",
       "      <td>2.122926e-13</td>\n",
       "      <td>0.073945</td>\n",
       "      <td>metabolite_100000010</td>\n",
       "      <td>Lachnospiraceae|Sellimonas</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>100 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                           r             p  \\\n",
       "(metabolite_100002253, Family_XIII|Family_XIII_...  0.508087  5.617983e-60   \n",
       "(metabolite_100002488, Ruminococcaceae|Ruminoco... -0.476256  7.957527e-54   \n",
       "(metabolite_100002488, Ruminococcaceae|Flavonif...  0.475678  1.108531e-53   \n",
       "(metabolite_100002488, Lachnospiraceae|Lachnocl...  0.468988  4.900222e-52   \n",
       "(metabolite_100000010, Family_XIII|Family_XIII_...  0.442081  4.059205e-43   \n",
       "...                                                      ...           ...   \n",
       "(metabolite_100020536, Lachnospiraceae|Lachnocl... -0.272563  9.977384e-17   \n",
       "(metabolite_100001417, Lachnospiraceae|Lachnosp... -0.265274  1.484930e-16   \n",
       "(metabolite_100001083, Erysipelotrichaceae|Erys... -0.263484  1.863074e-16   \n",
       "(metabolite_100000014, Ruminococcaceae|Ruminoco...  0.260556  2.574030e-16   \n",
       "(metabolite_100000010, Lachnospiraceae|Sellimonas) -0.271929  2.783801e-16   \n",
       "\n",
       "                                                        n             q  \\\n",
       "(metabolite_100002253, Family_XIII|Family_XIII_...  896.0  4.284274e-55   \n",
       "(metabolite_100002488, Ruminococcaceae|Ruminoco...  930.0  2.817885e-49   \n",
       "(metabolite_100002488, Ruminococcaceae|Flavonif...  930.0  2.817885e-49   \n",
       "(metabolite_100002488, Lachnospiraceae|Lachnocl...  930.0  9.342274e-48   \n",
       "(metabolite_100000010, Family_XIII|Family_XIII_...  874.0  6.191099e-39   \n",
       "...                                                   ...           ...   \n",
       "(metabolite_100020536, Lachnospiraceae|Lachnocl...  896.0  7.925785e-14   \n",
       "(metabolite_100001417, Lachnospiraceae|Lachnosp...  937.0  1.167431e-13   \n",
       "(metabolite_100001083, Erysipelotrichaceae|Erys...  944.0  1.449776e-13   \n",
       "(metabolite_100000014, Ruminococcaceae|Ruminoco...  957.0  1.982783e-13   \n",
       "(metabolite_100000010, Lachnospiraceae|Sellimonas)  874.0  2.122926e-13   \n",
       "\n",
       "                                                          r2  \\\n",
       "(metabolite_100002253, Family_XIII|Family_XIII_...  0.258152   \n",
       "(metabolite_100002488, Ruminococcaceae|Ruminoco...  0.226820   \n",
       "(metabolite_100002488, Ruminococcaceae|Flavonif...  0.226270   \n",
       "(metabolite_100002488, Lachnospiraceae|Lachnocl...  0.219949   \n",
       "(metabolite_100000010, Family_XIII|Family_XIII_...  0.195436   \n",
       "...                                                      ...   \n",
       "(metabolite_100020536, Lachnospiraceae|Lachnocl...  0.074291   \n",
       "(metabolite_100001417, Lachnospiraceae|Lachnosp...  0.070370   \n",
       "(metabolite_100001083, Erysipelotrichaceae|Erys...  0.069424   \n",
       "(metabolite_100000014, Ruminococcaceae|Ruminoco...  0.067889   \n",
       "(metabolite_100000010, Lachnospiraceae|Sellimonas)  0.073945   \n",
       "\n",
       "                                                              metabolite  \\\n",
       "(metabolite_100002253, Family_XIII|Family_XIII_...  metabolite_100002253   \n",
       "(metabolite_100002488, Ruminococcaceae|Ruminoco...  metabolite_100002488   \n",
       "(metabolite_100002488, Ruminococcaceae|Flavonif...  metabolite_100002488   \n",
       "(metabolite_100002488, Lachnospiraceae|Lachnocl...  metabolite_100002488   \n",
       "(metabolite_100000010, Family_XIII|Family_XIII_...  metabolite_100000010   \n",
       "...                                                                  ...   \n",
       "(metabolite_100020536, Lachnospiraceae|Lachnocl...  metabolite_100020536   \n",
       "(metabolite_100001417, Lachnospiraceae|Lachnosp...  metabolite_100001417   \n",
       "(metabolite_100001083, Erysipelotrichaceae|Erys...  metabolite_100001083   \n",
       "(metabolite_100000014, Ruminococcaceae|Ruminoco...  metabolite_100000014   \n",
       "(metabolite_100000010, Lachnospiraceae|Sellimonas)  metabolite_100000010   \n",
       "\n",
       "                                                                                            taxon  \n",
       "(metabolite_100002253, Family_XIII|Family_XIII_...                Family_XIII|Family_XIII_UCG-001  \n",
       "(metabolite_100002488, Ruminococcaceae|Ruminoco...  Ruminococcaceae|Ruminococcaceae_NK4A214_group  \n",
       "(metabolite_100002488, Ruminococcaceae|Flavonif...                 Ruminococcaceae|Flavonifractor  \n",
       "(metabolite_100002488, Lachnospiraceae|Lachnocl...              Lachnospiraceae|Lachnoclostridium  \n",
       "(metabolite_100000010, Family_XIII|Family_XIII_...                Family_XIII|Family_XIII_UCG-001  \n",
       "...                                                                                           ...  \n",
       "(metabolite_100020536, Lachnospiraceae|Lachnocl...              Lachnospiraceae|Lachnoclostridium  \n",
       "(metabolite_100001417, Lachnospiraceae|Lachnosp...        Lachnospiraceae|Lachnospiraceae_UCG-004  \n",
       "(metabolite_100001083, Erysipelotrichaceae|Erys...     Erysipelotrichaceae|Erysipelatoclostridium  \n",
       "(metabolite_100000014, Ruminococcaceae|Ruminoco...        Ruminococcaceae|Ruminococcaceae_UCG-002  \n",
       "(metabolite_100000010, Lachnospiraceae|Sellimonas)                     Lachnospiraceae|Sellimonas  \n",
       "\n",
       "[100 rows x 7 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metrics.sort_values(by=\"p\").head(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "amateur-costs",
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics[metrics.q < 0.05].to_csv(\"data/sig_metabolite_taxon.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cardiac-director",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "508"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metrics[metrics.q < 0.05].metabolite.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indirect-status",
   "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
}