GLM_model_on_intersample_concentration_ratio_for_targeted_compounds.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "004e916c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paths import *\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import statsmodels.api as sm\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from tqdm import tqdm\n",
    "import joblib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import seaborn as sns\n",
    "import multiprocessing as mp\n",
    "import warnings\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "893ddf60",
   "metadata": {},
   "outputs": [],
   "source": [
    "NCORES = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "32dcdbe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "warnings.filterwarnings(\"ignore\", message=\"`distplot` is a deprecated function\")\n",
    "warnings.filterwarnings(\"ignore\", message=\"DeprecationWarning\")\n",
    "\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "deae1e40",
   "metadata": {},
   "outputs": [],
   "source": [
    "def glm_model(data_, analyte_):\n",
    "    data_analyte = data_[data_[\"Analyte\"] == analyte_]\n",
    "    le = LabelEncoder()\n",
    "    data_analyte.loc[:, 'Disease_label'] = le.fit_transform(data_analyte['Group Name'])\n",
    "    data_analyte.loc[:, 'Gender_label'] = le.fit_transform(data_analyte['Gender'])\n",
    "    data_analyte.loc[:, 'site'] = data_analyte.loc[:, 'Client Sample ID'].str[:3]\n",
    "    data_analyte = data_analyte[[\"Result\", \"Disease_label\", \"Age\", \"BMI\", \"Gender_label\", \"site\", \"house\"]]\n",
    "    X = sm.add_constant(data_analyte[['Age', 'BMI', 'Gender_label', 'Disease_label']])\n",
    "    site_dummies = pd.get_dummies(data_analyte['site'], prefix='site', drop_first=True)\n",
    "    house_dummies = pd.get_dummies(data_analyte['house'], prefix='house', drop_first=True)\n",
    "    X = pd.concat([X, site_dummies, house_dummies], axis=1)\n",
    "    model = sm.GLM(data_analyte['Result'], X, family=sm.families.Gamma(link=sm.families.links.log()))\n",
    "#     model = sm.GLM(data_analyte['Result'], X, family=sm.families.Gaussian(link=sm.families.links.identity()))\n",
    "    mdf = model.fit()\n",
    "    out_dict = {}\n",
    "    out_dict[\"analyte\"] = analyte_\n",
    "    out_dict[\"sample\"] = \"serum_feces_ratio\"\n",
    "    out_dict[\"patient_sample_count\"] = data_analyte.shape[0]\n",
    "    out_dict[\"model\"] = mdf\n",
    "    return out_dict \n",
    "\n",
    "def get_model_parallel(data_, analytes_):\n",
    "    data_list = [data_]*len(analytes_)\n",
    "    analytes_list = list(analytes_)\n",
    "    arg_list = list(zip(data_list, analytes_list))\n",
    "    p = mp.Pool(NCORES)\n",
    "    model_list = p.starmap(glm_model, arg_list)\n",
    "    p.close()\n",
    "    p.join()\n",
    "    return model_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "33f618dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "filename = SHORT_CHAIN_FATTY_ACID_DATA_FILENAME\n",
    "mapping_filename = \"short_chain_fatty_acid_spoke_map.csv\"\n",
    "file_path = os.path.join(DATA_ROOT_PATH, filename)\n",
    "mapping_filepath = os.path.join(OUTPUT_PATH, mapping_filename)\n",
    "data = pd.read_excel(file_path, engine='openpyxl')\n",
    "data_feces = data[data[\"Client Matrix\"]==\"feces\"]\n",
    "data_serum = data[data[\"Client Matrix\"]==\"serum\"]\n",
    "mapping_data = pd.read_csv(mapping_filepath)\n",
    "analytes = mapping_data[\"name\"].unique()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "4b661417",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/frame.py:4174: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  errors=errors,\n"
     ]
    }
   ],
   "source": [
    "data_feces.loc[:, 'house'] = (data_feces['Client Sample ID'].str[:3] + data_feces['Client Sample ID'].str[-4:])\n",
    "data_feces_exclude_outlier_threshold_column = data_feces.drop(\"Analysis Comment\", axis=1)\n",
    "house_to_exclude = data_feces_exclude_outlier_threshold_column[data_feces_exclude_outlier_threshold_column.isna().any(axis=1)].house.values\n",
    "data_feces_nan_removed = data_feces[~data_feces[\"house\"].isin(house_to_exclude)]\n",
    "data_feces_outlier_removed = data_feces_nan_removed[data_feces_nan_removed[\"Analysis Comment\"].isna()]\n",
    "data_feces_outlier_removed.drop(\"Analysis Comment\", axis=1, inplace=True)\n",
    "group_counts = data_feces_outlier_removed.groupby(['Analyte', 'house']).size()\n",
    "valid_index = group_counts[group_counts == 2].reset_index().set_index(['Analyte', 'house']).index\n",
    "data_feces_outlier_removed = data_feces_outlier_removed.set_index(['Analyte', 'house']).loc[valid_index].reset_index()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "e4c2ec3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/frame.py:4174: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  errors=errors,\n"
     ]
    }
   ],
   "source": [
    "data_serum.loc[:, 'house'] = (data_serum['Client Sample ID'].str[:3] + data_serum['Client Sample ID'].str[-4:])\n",
    "data_serum_exclude_outlier_threshold_column = data_serum.drop(\"Analysis Comment\", axis=1)\n",
    "house_to_exclude = data_serum_exclude_outlier_threshold_column[data_serum_exclude_outlier_threshold_column.isna().any(axis=1)].house.values\n",
    "data_serum_nan_removed = data_serum[~data_serum[\"house\"].isin(house_to_exclude)]\n",
    "data_serum_outlier_removed = data_serum_nan_removed[data_serum_nan_removed[\"Analysis Comment\"].isna()]\n",
    "data_serum_outlier_removed.drop(\"Analysis Comment\", axis=1, inplace=True)\n",
    "group_counts = data_serum_outlier_removed.groupby(['Analyte', 'house']).size()\n",
    "valid_index = group_counts[group_counts == 2].reset_index().set_index(['Analyte', 'house']).index\n",
    "data_serum_outlier_removed = data_serum_outlier_removed.set_index(['Analyte', 'house']).loc[valid_index].reset_index()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "003ec000",
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_list = list(data_serum_outlier_removed.columns)\n",
    "columns_to_remove = [\"Unique Tube Label ID\", \"Dilution Factor\", \"Unit\", \"Race / Ethnicity\", \"Sample Description\", \"Client Matrix\"]\n",
    "\n",
    "columns_to_merge = [\"Client Sample ID\", \"Analyte\", \"Group Name\"]\n",
    "columns_list_x = list(map(lambda x:x+\"_x\", columns_list))\n",
    "for index, item in enumerate(columns_list):\n",
    "    if item in columns_to_remove:\n",
    "        columns_list_x.remove(item+\"_x\")\n",
    "for item in columns_to_merge:\n",
    "    columns_list_x.remove(item+\"_x\")\n",
    "    columns_list_x.append(item)\n",
    "    \n",
    "data_serum_outlier_removed.drop(columns_to_remove, axis=1, inplace=True)\n",
    "data_feces_outlier_removed.drop(columns_to_remove, axis=1, inplace=True)\n",
    "\n",
    "data_serum_feces_outlier_removed = pd.merge(data_serum_outlier_removed, data_feces_outlier_removed, on=columns_to_merge)\n",
    "data_serum_feces_outlier_removed.loc[:, \"Result\"] = np.divide(data_serum_feces_outlier_removed.Result_x, data_serum_feces_outlier_removed.Result_y)\n",
    "data_serum_feces_outlier_removed = data_serum_feces_outlier_removed[[\"Result\"] + columns_list_x]\n",
    "data_serum_feces_outlier_removed.drop(\"Result_x\", axis=1, inplace=True)\n",
    "\n",
    "for item in columns_list:\n",
    "    if item+\"_x\" in list(data_serum_feces_outlier_removed.columns):\n",
    "        data_serum_feces_outlier_removed.rename(columns={item+\"_x\":item}, inplace=True)\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "adc65d16",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1596: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[key] = _infer_fill_value(value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/indexing.py:1743: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  isetter(ilocs[0], value)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/generalized_linear_model.py:798: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  return np.sum(resid / self.family.variance(mu)) / self.df_resid\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 126 ms, sys: 133 ms, total: 259 ms\n",
      "Wall time: 9.11 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['/pool0/home/karthik/iMSMS_metabolomics/result/GLM_shortChain_fattyAcid_without_outlier_sample_serum_feces_ratio.joblib']"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "model_without_outlier_list = get_model_parallel(data_serum_feces_outlier_removed, analytes)\n",
    "joblib.dump(model_without_outlier_list, os.path.join(OUTPUT_PATH, \"GLM_shortChain_fattyAcid_without_outlier_sample_serum_feces_ratio.joblib\"))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13ae90cd",
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}