GLM_model_on_intersample_concentration_ratio_for_global_compounds.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "d7b32e8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paths import *\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "import statsmodels.api as sm\n",
    "import multiprocessing as mp\n",
    "import sys\n",
    "import time\n",
    "import joblib\n",
    "from tqdm import tqdm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "2c9c8f8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "data_type = \"without_outlier\"\n",
    "NCORES = 12\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "0941d5dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def glm_model(data_, analyte_, remove_outlier):\n",
    "    data_analyte = data_[selected_metadata_columns+[analyte_, \"site\", \"house\"]]\n",
    "    if remove_outlier:\n",
    "        data_analyte = remove_outliers(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 = data_analyte[[analyte_, \"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",
    "    try:\n",
    "        model = sm.GLM(data_analyte[analyte_], X, family=sm.families.Gamma(link=sm.families.links.log()))\n",
    "        mdf = model.fit()\n",
    "    except:\n",
    "        mdf = None\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_, remove_outlier_flag):    \n",
    "    data_list = [data_]*len(analytes_)\n",
    "    analytes_list = list(analytes_)\n",
    "    remove_outlier_flag_list = [remove_outlier_flag]*len(analytes_)\n",
    "    arg_list = list(zip(data_list, analytes_list, remove_outlier_flag_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\n",
    "\n",
    "def remove_outliers(data_, analyte):\n",
    "    analyte_concentration_values = data_[analyte]\n",
    "    # Calculate quartiles and IQR\n",
    "    q1 = analyte_concentration_values.quantile(0.25)\n",
    "    q3 = analyte_concentration_values.quantile(0.75)\n",
    "    iqr = q3 - q1\n",
    "    # Define Tukey's fences\n",
    "    k = 3\n",
    "    lower_fence = q1 - k * iqr\n",
    "    upper_fence = q3 + k * iqr\n",
    "    # Filter non-outliers\n",
    "    non_outliers = data_[(data_[analyte] >= lower_fence) & (data_[analyte] <= upper_fence)]\n",
    "    house_counts = non_outliers.groupby(\"house\").size()\n",
    "    # Filter to keep only the rows with paired house IDs\n",
    "    paired_houses = house_counts[house_counts == 2].index\n",
    "    non_outliers = non_outliers[non_outliers[\"house\"].isin(paired_houses)]\n",
    "    return non_outliers\n",
    "\n",
    "def load_data(sample):\n",
    "    if sample == \"serum\":\n",
    "        filename = GLOBAL_SERUM_DATA_FILENAME\n",
    "    else:\n",
    "        filename = GLOBAL_STOOL_DATA_FILENAME\n",
    "    file_path = os.path.join(DATA_ROOT_PATH, filename)\n",
    "    sheet_name = [\"Chemical Annotation\", \"Sample Meta Data\", \"Log Transformed Data\"]\n",
    "    analyte_metadata = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[0])\n",
    "    patient_metadata = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[1])\n",
    "    data = pd.read_excel(file_path, engine='openpyxl', sheet_name=sheet_name[2])\n",
    "    return analyte_metadata, patient_metadata, data\n",
    "\n",
    "\n",
    "def reverse_log_transform(data_, analyte_columns_selected_):\n",
    "    data_with_patient_metadata_reverse_log_transformed = data_.copy()\n",
    "    for analyte in analyte_columns_selected_:\n",
    "        analyte_concentration = data_with_patient_metadata_reverse_log_transformed[analyte].values\n",
    "        data_with_patient_metadata_reverse_log_transformed.loc[:, analyte] = np.exp(analyte_concentration)\n",
    "    return data_with_patient_metadata_reverse_log_transformed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "18833513",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_metadata_columns = [\"PARENT_SAMPLE_NAME\", \"CLIENT_IDENTIFIER\", \"GROUP_NAME\", \"AGE\", \"BMI\", \"GENDER\", \"CLIENT_SAMPLE_ID\", \"CLIENT_MATRIX\", \"TREATMENT\"]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "6a81d670",
   "metadata": {},
   "outputs": [],
   "source": [
    "analyte_metadata_feces, patient_metadata_feces, data_feces = load_data(\"feces\")\n",
    "analyte_metadata_serum, patient_metadata_serum, data_serum = load_data(\"serum\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "eb78ee62",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_data(sample):\n",
    "    analyte_metadata, patient_metadata, data = load_data(sample)\n",
    "    global_metabolomics_compound_spoke_map = pd.read_csv(os.path.join(OUTPUT_PATH, \"global_metabolomics_compound_spoke_map.csv\"))\n",
    "    analyte_columns = list(data.columns)\n",
    "    analyte_columns.remove(\"PARENT_SAMPLE_NAME\")\n",
    "    analyte_columns_selected = global_metabolomics_compound_spoke_map[global_metabolomics_compound_spoke_map.CHEM_ID.isin(analyte_columns)][\"CHEM_ID\"].unique()\n",
    "    data_with_analyte_columns_selected = data[[\"PARENT_SAMPLE_NAME\"]+list(analyte_columns_selected)]    \n",
    "    patient_metadata_selected_columns = patient_metadata[selected_metadata_columns]\n",
    "    patient_metadata_selected_columns.loc[:, 'house'] = (patient_metadata_selected_columns['CLIENT_SAMPLE_ID'].str[:3] + patient_metadata_selected_columns['CLIENT_SAMPLE_ID'].str[-4:])\n",
    "    patient_metadata_selected_columns.loc[:, 'site'] = patient_metadata_selected_columns.loc[:, 'CLIENT_SAMPLE_ID'].str[:3]\n",
    "    le = LabelEncoder()\n",
    "    patient_metadata_selected_columns.loc[:, 'Disease_label'] = le.fit_transform(patient_metadata_selected_columns['GROUP_NAME'])\n",
    "    patient_metadata_selected_columns.loc[:, 'Gender_label'] = le.fit_transform(patient_metadata_selected_columns['GENDER'])\n",
    "    data_with_patient_metadata = pd.merge(data_with_analyte_columns_selected, patient_metadata_selected_columns, on=\"PARENT_SAMPLE_NAME\")\n",
    "    # Reverse log transforming the data, so that gamma family can be applied in GLM model\n",
    "    data_with_patient_metadata_reverse_log_transformed = reverse_log_transform(data_with_patient_metadata, analyte_columns_selected)\n",
    "    #Removing patients (and their partners) with missing data\n",
    "    house_to_exclude = data_with_patient_metadata_reverse_log_transformed[data_with_patient_metadata_reverse_log_transformed.isna().any(axis=1)].house.values\n",
    "    data_with_patient_metadata_reverse_log_transformed_nan_removed = data_with_patient_metadata_reverse_log_transformed[~data_with_patient_metadata_reverse_log_transformed[\"house\"].isin(house_to_exclude)]\n",
    "    data_to_go = data_with_patient_metadata_reverse_log_transformed_nan_removed\n",
    "    if data_type == \"with_outlier\":\n",
    "        data_out = data_to_go  \n",
    "    elif data_type == \"with_outlier_treated\":\n",
    "        data_to_go_treated_ = data_to_go[data_to_go.TREATMENT == \"Treated\"]\n",
    "        house_to_include = data_to_go_treated_[\"house\"].unique()\n",
    "        data_to_go_treated = data_to_go[data_to_go.house.isin(house_to_include)]\n",
    "        data_out = data_to_go_treated\n",
    "    elif data_type == \"with_outlier_not_treated\":\n",
    "        data_to_go_not_treated_ = data_to_go[data_to_go.TREATMENT == \"Off\"]\n",
    "        house_to_include = data_to_go_not_treated_[\"house\"].unique()\n",
    "        data_to_go_not_treated = data_to_go[data_to_go.house.isin(house_to_include)]\n",
    "        data_out = data_to_go_not_treated\n",
    "    elif data_type == \"without_outlier\":\n",
    "        data_out = data_to_go\n",
    "    elif data_type == \"without_outlier_treated\":\n",
    "        data_to_go_treated_ = data_to_go[data_to_go.TREATMENT == \"Treated\"]\n",
    "        house_to_include = data_to_go_treated_[\"house\"].unique()\n",
    "        data_to_go_treated = data_to_go[data_to_go.house.isin(house_to_include)]\n",
    "        data_out = data_to_go_treated\n",
    "    elif data_type == \"without_outlier_not_treated\":\n",
    "        data_to_go_not_treated_ = data_to_go[data_to_go.TREATMENT == \"Off\"]\n",
    "        house_to_include = data_to_go_not_treated_[\"house\"].unique()\n",
    "        data_to_go_not_treated = data_to_go[data_to_go.house.isin(house_to_include)]\n",
    "        data_out = data_to_go_not_treated\n",
    "    return data_out\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "360b8900",
   "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: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",
      "100%|██████████| 1206/1206 [00:05<00:00, 213.63it/s]\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/pandas/core/frame.py:4308: 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 = get_data(\"feces\")\n",
    "data_serum = get_data(\"serum\")\n",
    "\n",
    "non_compound_columns = [\"PARENT_SAMPLE_NAME\", \"AGE\", \"BMI\", \"GENDER\", \"CLIENT_SAMPLE_ID\", \"CLIENT_MATRIX\", \"TREATMENT\", \"house\", \"site\", \"Disease_label\", \"Gender_label\", \"CLIENT_IDENTIFIER\", \"GROUP_NAME\"]\n",
    "non_compound_columns_x = list(map(lambda x:x+\"_x\", non_compound_columns))\n",
    "non_compound_columns_y = list(map(lambda x:x+\"_y\", non_compound_columns))\n",
    "non_compound_columns_x.pop(4)\n",
    "non_compound_columns_x.append(\"CLIENT_SAMPLE_ID\")\n",
    "serum_compound_columns = list(data_serum.columns[~data_serum.columns.isin(non_compound_columns)])\n",
    "\n",
    "data_serum_feces = pd.merge(data_serum, data_feces, on=\"CLIENT_SAMPLE_ID\")\n",
    "new_compound_columns = []\n",
    "for serum_compound_column in tqdm(serum_compound_columns):    \n",
    "    try:\n",
    "        new_column_name = serum_compound_column\n",
    "        data_serum_feces.loc[:, new_column_name] = np.divide(data_serum_feces[str(serum_compound_column)+\"_x\"], data_serum_feces[str(serum_compound_column)+\"_y\"])\n",
    "        new_compound_columns.append(new_column_name)\n",
    "    except KeyError:\n",
    "        continue\n",
    "    \n",
    "data_serum_feces_to_go = data_serum_feces[new_compound_columns+non_compound_columns_x]\n",
    "for item in non_compound_columns:\n",
    "    if item != \"CLIENT_SAMPLE_ID\":\n",
    "        data_serum_feces_to_go.rename(columns={item+\"_x\":item}, inplace=True)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "5835a19d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n",
      "/pool0/home/karthik/iMSMS_metabolomics/venv/lib/python3.6/site-packages/statsmodels/genmod/families/family.py:179: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.sum(resid_dev * freq_weights * var_weights / scale)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.56 s, sys: 8.78 s, total: 12.3 s\n",
      "Wall time: 22min 43s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "global_metabolomics_compound_spoke_map = pd.read_csv(os.path.join(OUTPUT_PATH, \"global_metabolomics_compound_spoke_map.csv\"))\n",
    "analyte_columns_feces = list(data_feces.columns)\n",
    "analyte_columns_feces.remove(\"PARENT_SAMPLE_NAME\")\n",
    "analyte_columns_feces_selected = global_metabolomics_compound_spoke_map[global_metabolomics_compound_spoke_map.CHEM_ID.isin(analyte_columns_feces)][\"CHEM_ID\"].unique()\n",
    "\n",
    "analyte_columns_serum = list(data_serum.columns)\n",
    "analyte_columns_serum.remove(\"PARENT_SAMPLE_NAME\")\n",
    "analyte_columns_serum_selected = global_metabolomics_compound_spoke_map[global_metabolomics_compound_spoke_map.CHEM_ID.isin(analyte_columns_serum)][\"CHEM_ID\"].unique()\n",
    "\n",
    "analyte_columns_selected = list(set(analyte_columns_serum_selected).intersection(set(analyte_columns_feces_selected)))\n",
    "\n",
    "\n",
    "model_list = get_model_parallel(data_serum_feces_to_go, analyte_columns_selected, True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "aa36600d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/pool0/home/karthik/iMSMS_metabolomics/result/GLM_global_compounds_without_outlier_sample_serum_feces_ratio.joblib']"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_filename = \"GLM_global_compounds_without_outlier_sample_serum_feces_ratio.joblib\"\n",
    "joblib.dump(model_list, os.path.join(OUTPUT_PATH, model_filename))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3ab6e0f",
   "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
}