GRPM_03_dataset_survey.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GRPM MeSH Screening \n",
    "\n",
    "This notebook is engineered to screen the previously retrieved genetic polymorphism data using selected MeSH terms. It works with MeSH sets that are used as hooks to retrieve subsets of genes and polymorphisms from the \"GRPM ds\" dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Only for Google Colab\n",
    "import os\n",
    "import sys\n",
    "import subprocess\n",
    "\n",
    "# @markdown Run in Colab virtual machine by default\n",
    "\n",
    "# @markdown to run in google drive set:\n",
    "import_mydrive = False #@param {type:\"boolean\"}\n",
    "\n",
    "\n",
    "if 'google.colab' in sys.modules:\n",
    "    subprocess.run([\"pip\", \"install\", \"nbib\"])\n",
    "    subprocess.run([\"pip\", \"install\", \"biopython\"])\n",
    "\n",
    "    if import_mydrive:\n",
    "        from google.colab import drive\n",
    "        drive.mount('/content/drive')\n",
    "\n",
    "        if not os.path.exists('/content/drive/MyDrive/grpm_system/'):\n",
    "            subprocess.run(['mkdir', '/content/drive/MyDrive/grpm_system/'])\n",
    "        subprocess.run(['cd', '/content/drive/MyDrive/grpm_system/'])\n",
    "    else:\n",
    "        if not os.path.exists('/content/grpm_system/'):\n",
    "            subprocess.run(['mkdir', '/content/grpm_system/'])\n",
    "        subprocess.run(['cd', '/content/grpm_system/'])\n",
    "\n",
    "current_directory = os.getcwd()\n",
    "print(\"Current working directory:\", current_directory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T13:51:23.567025800Z",
     "start_time": "2024-04-13T13:51:23.262849600Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import io\n",
    "import glob\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import requests\n",
    "import zipfile\n",
    "\n",
    "def simple_bool(message):\n",
    "    choose = input(message+\" (y/n): \").lower()\n",
    "    your_bool = choose in [\"y\", \"yes\",\"yea\",\"sure\"]\n",
    "    return your_bool\n",
    "\n",
    "def get_and_extract(file, dir = os.getcwd()):\n",
    "    url = \"https://zenodo.org/record/8205724/files/\"+file+\".zip?download=1\"\n",
    "    zip_file_name = file+\".zip\"\n",
    "    extracted_folder_name = dir\n",
    "\n",
    "    # Download the ZIP file\n",
    "    response = requests.get(url)\n",
    "\n",
    "    if response.status_code == 200:\n",
    "        # Extract the ZIP contents\n",
    "        with io.BytesIO(response.content) as zip_buffer:\n",
    "            with zipfile.ZipFile(zip_buffer, 'r') as zip_ref:\n",
    "                zip_ref.extractall(extracted_folder_name)\n",
    "        print(f\"ZIP file '{zip_file_name}' extracted to '{extracted_folder_name}' successfully.\")\n",
    "    else:\n",
    "        print(\"Failed to download the ZIP file.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get requirements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Get GRPM Dataset from Zenodo Repository\n",
    "#https://zenodo.org/record/8205724  DOI: 10.5281/zenodo.8205724\n",
    "\n",
    "if simple_bool('Download pre-made GRPM-Dataset from Zenodo? (6.5 minutes in Colab)'):\n",
    "    timea = datetime.now()\n",
    "    get_and_extract('grpm_dataset')\n",
    "    print('Download and extraction time ',datetime.now()-timea)\n",
    "\n",
    "if simple_bool('Download pre-made ref-mesh-archive from Zenodo?'):\n",
    "    timea = datetime.now()\n",
    "    get_and_extract('ref-mesh-archive')\n",
    "    print('Download and extraction time ',datetime.now()-timea)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import GRPM dataset (required)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T13:52:28.037954400Z",
     "start_time": "2024-04-13T13:51:23.557029Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "importing GRPM Dataset...\n",
      "time load: 0:01:02.019286\n"
     ]
    }
   ],
   "source": [
    "#Load GRPM db Report-----------------------------------------\n",
    "\n",
    "# choose database:\n",
    "db_tag      = 'pcg'\n",
    "# 'pcg'    = protein coding genes = grpm_db\n",
    "# 'rna'    = rna genes            = grpm_db_rna\n",
    "# 'pseudo' = pseudogenes          = grpm_db_pseudo\n",
    "\n",
    "db_name = 'grpm_db_'+ db_tag\n",
    "db_path = 'grpm_dataset/'+db_name\n",
    "\n",
    "print('importing GRPM Dataset...')\n",
    "#get gene list from grpm report\n",
    "GRPM_report = pd.read_csv(db_path+'/GRPM_report.csv',index_col=0).transpose().reset_index().rename(columns={'index':'gene'})\n",
    "grpm_genes_list = GRPM_report.gene.to_list()\n",
    "\n",
    "#Import grpm data back-------------------------------------------\n",
    "time_load_1 = datetime.now()\n",
    "\n",
    "columns = ['gene', 'rsid', 'pmids', 'mesh']\n",
    "grpm_dataset = pd.read_csv(db_path+'/grpm_table_output.csv', usecols=columns)\n",
    "\n",
    "grpm_dataset['pmids'] = grpm_dataset['pmids'].astype(str) #convert pmid type in str\n",
    "time_load_2 = datetime.now()\n",
    "print('time load:',time_load_2-time_load_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subset GRPM Dataset (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subset_grpm = simple_bool('Do you want to use a custom gene list to subset GRPM Dataset?')\n",
    "if subset_grpm:\n",
    "    # import your custom gene list (.csv)\n",
    "    file_csv = []\n",
    "    for file in os.listdir():\n",
    "        if file.endswith(\".csv\") or file.endswith(\".tsv\"):\n",
    "            file_csv.append(file)\n",
    "\n",
    "    filenum = input('import your custom gene list (.csv)\\nselect file index: \\n'+str(pd.Series(file_csv)))\n",
    "\n",
    "    time1 = datetime.now()\n",
    "    subset_genes = pd.read_csv(file_csv[int(filenum)])\n",
    "    subset_genes = subset_genes[subset_genes.columns[int(input('select column index:\\n'+ str(pd.Series(subset_genes.columns))))]].drop_duplicates().str.replace(' ','')\n",
    "    subset_genes.to_list()\n",
    "\n",
    "    # subsetting GRPM_report and grpm_dataset\n",
    "    GRPM_report_subset = GRPM_report[GRPM_report['gene'].isin(subset_genes)]\n",
    "    grpm_subset = grpm_dataset[grpm_dataset['gene'].isin(subset_genes)]\n",
    "    print(\"You're using a GRPM Dataset partition\\ntime subsetting:\",time_load_2-time_load_1)\n",
    "    display(GRPM_report_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define context\n",
    "    - gene list\n",
    "    - survey directory\n",
    "    - ref-mesh list"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Check avalable ref-MeSH lists\n",
    "## Set directory/import data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available reference mesh lists:\n",
      "\n",
      " 0             Physical Endurance\n",
      "1              Nutrition Therapy\n",
      "5               Nutrition Policy\n",
      "6                 Motor Activity\n",
      "7          Adipose Tissue, Brown\n",
      "                 ...            \n",
      "623    Body Weights and Measures\n",
      "627                  Body Weight\n",
      "628      Body Weight Maintenance\n",
      "629                   Overweight\n",
      "630                      Obesity\n",
      "Name: mesh, Length: 243, dtype: object\n"
     ]
    }
   ],
   "source": [
    "#Check avalable refs:\n",
    "ref_path = \"ref-mesh-archive/\"  # Replace with the actual ref mesh path\n",
    "\n",
    "#---------------------------------\n",
    "#use random mesh list?\n",
    "random_mesh = False\n",
    "if random_mesh:\n",
    "    ref_path = \"ref-mesh-archive/random_lists/\"\n",
    "#---------------------------------\n",
    "\n",
    "# Create a file path pattern to match CSV files\n",
    "file_pattern = os.path.join(ref_path, \"*.csv\")\n",
    "\n",
    "# Use glob to get a list of file paths matching the pattern\n",
    "csv_files = glob.glob(file_pattern)\n",
    "\n",
    "csv_files_name = []\n",
    "# Print the list of CSV files\n",
    "for file in csv_files:\n",
    "    file_name = os.path.basename(file)\n",
    "    csv_files_name.append(file_name)\n",
    "\n",
    "pd.set_option('display.max_rows', 100)\n",
    "print('Available reference mesh lists:')\n",
    "csv_files_df = pd.Series(csv_files_name)\n",
    "\n",
    "csv_file_tag = pd.DataFrame()\n",
    "if not random_mesh:\n",
    "    csv_file_tag = csv_files_df.str.extract(r'ref_mesh_(.*)\\.csv', expand=False).dropna().reset_index(drop=True)\n",
    "else:\n",
    "    csv_file_tag = csv_files_df.str.extract(r'(.*)\\.csv', expand=False).dropna().reset_index(drop=True)\n",
    "\n",
    "\n",
    "#------------------------------------------------------\n",
    "# define directory folder path:\n",
    "survey_path = 'grpm_surveys/' # keep default to use root path\n",
    "\n",
    "# choose ref_mesh.csv tab:\n",
    "topic_tag   = csv_file_tag[int(input('\\Select index from available ref-mesh list:\\n'+str(csv_file_tag)))]\n",
    "add         = ''    # additional survey directory tag\n",
    "#------------------------------------------------------\n",
    "\n",
    "\n",
    "# (1) Create survey directory:\n",
    "survey_path = survey_path+'grpm_random/' if random_mesh else survey_path\n",
    "directory = survey_path + 'grpm_survey_' + db_tag + '_' + topic_tag + add\n",
    "if not os.path.exists(directory):\n",
    "    os.makedirs(directory)\n",
    "\n",
    "\n",
    "# (2) Import Mesh-reference list:\n",
    "ref_filename = \"ref_mesh_\" + topic_tag + \".csv\" if not random_mesh else topic_tag + \".csv\"\n",
    "ref = pd.read_csv(ref_path + ref_filename, index_col=0)\n",
    "\n",
    "if 'mesh' not in ref.columns:\n",
    "    ref = ref.rename(columns={'Preferred Label': 'mesh'})\n",
    "\n",
    "ref_mesh_n = ref.mesh.nunique()\n",
    "ref_mesh_list = ref['mesh'].drop_duplicates()\n",
    "\n",
    "#----------------------------------------------------------\n",
    "print('\\n', ref_mesh_list)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T13:53:39.450710900Z",
     "start_time": "2024-04-13T13:53:29.030415Z"
    }
   },
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Get GRPM dataset Metrics [required]\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: total: 0 ns\n",
      "Wall time: 0 ns\n",
      "0:00:03.064639\n",
      "0:01:09.662249\n"
     ]
    },
    {
     "data": {
      "text/plain": "         rsid                          pmids                        mesh  \\\n        count unique          top freq count unique       top freq count   \ngene                                                                       \nA1BG      118      4     rs893184   68   118      7  23690342   40   118   \nA1CF      362      7   rs10821905  185   362     23  29437585   31   362   \nA2M      1160     21        rs669  676  1160     57  32747830   92  1160   \nA2ML1    1065     32  rs863224951  119  1065     13  31009165  575  1065   \nA3GALT2    19      1  rs376200069   19    19      1  28506304   19    19   \n\n                                    pmidmesh         \\\n        unique             top freq    count unique   \ngene                                                  \nA1BG        75    Hypertension    8      118     94   \nA1CF       137          Humans   22      362    318   \nA2M        365          Humans   61     1160    887   \nA2ML1      117    Otitis Media   51     1065    227   \nA3GALT2     11  Adenocarcinoma    3       19     11   \n\n                                             \n                                   top freq  \ngene                                         \nA1BG                     32279138HIV-1    3  \nA1CF      30529582Colorectal Neoplasms    4  \nA2M      24011543Amyloid beta-Peptides    8  \nA2ML1               31009165Adolescent   23  \nA3GALT2         28506304Adenocarcinoma    3  ",
      "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 tr th {\n        text-align: left;\n    }\n\n    .dataframe thead tr:last-of-type th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr>\n      <th></th>\n      <th colspan=\"4\" halign=\"left\">rsid</th>\n      <th colspan=\"4\" halign=\"left\">pmids</th>\n      <th colspan=\"4\" halign=\"left\">mesh</th>\n      <th colspan=\"4\" halign=\"left\">pmidmesh</th>\n    </tr>\n    <tr>\n      <th></th>\n      <th>count</th>\n      <th>unique</th>\n      <th>top</th>\n      <th>freq</th>\n      <th>count</th>\n      <th>unique</th>\n      <th>top</th>\n      <th>freq</th>\n      <th>count</th>\n      <th>unique</th>\n      <th>top</th>\n      <th>freq</th>\n      <th>count</th>\n      <th>unique</th>\n      <th>top</th>\n      <th>freq</th>\n    </tr>\n    <tr>\n      <th>gene</th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>A1BG</th>\n      <td>118</td>\n      <td>4</td>\n      <td>rs893184</td>\n      <td>68</td>\n      <td>118</td>\n      <td>7</td>\n      <td>23690342</td>\n      <td>40</td>\n      <td>118</td>\n      <td>75</td>\n      <td>Hypertension</td>\n      <td>8</td>\n      <td>118</td>\n      <td>94</td>\n      <td>32279138HIV-1</td>\n      <td>3</td>\n    </tr>\n    <tr>\n      <th>A1CF</th>\n      <td>362</td>\n      <td>7</td>\n      <td>rs10821905</td>\n      <td>185</td>\n      <td>362</td>\n      <td>23</td>\n      <td>29437585</td>\n      <td>31</td>\n      <td>362</td>\n      <td>137</td>\n      <td>Humans</td>\n      <td>22</td>\n      <td>362</td>\n      <td>318</td>\n      <td>30529582Colorectal Neoplasms</td>\n      <td>4</td>\n    </tr>\n    <tr>\n      <th>A2M</th>\n      <td>1160</td>\n      <td>21</td>\n      <td>rs669</td>\n      <td>676</td>\n      <td>1160</td>\n      <td>57</td>\n      <td>32747830</td>\n      <td>92</td>\n      <td>1160</td>\n      <td>365</td>\n      <td>Humans</td>\n      <td>61</td>\n      <td>1160</td>\n      <td>887</td>\n      <td>24011543Amyloid beta-Peptides</td>\n      <td>8</td>\n    </tr>\n    <tr>\n      <th>A2ML1</th>\n      <td>1065</td>\n      <td>32</td>\n      <td>rs863224951</td>\n      <td>119</td>\n      <td>1065</td>\n      <td>13</td>\n      <td>31009165</td>\n      <td>575</td>\n      <td>1065</td>\n      <td>117</td>\n      <td>Otitis Media</td>\n      <td>51</td>\n      <td>1065</td>\n      <td>227</td>\n      <td>31009165Adolescent</td>\n      <td>23</td>\n    </tr>\n    <tr>\n      <th>A3GALT2</th>\n      <td>19</td>\n      <td>1</td>\n      <td>rs376200069</td>\n      <td>19</td>\n      <td>19</td>\n      <td>1</td>\n      <td>28506304</td>\n      <td>19</td>\n      <td>19</td>\n      <td>11</td>\n      <td>Adenocarcinoma</td>\n      <td>3</td>\n      <td>19</td>\n      <td>11</td>\n      <td>28506304Adenocarcinoma</td>\n      <td>3</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time\n",
    "# GET GRPM whole dataset Metrics\n",
    "time_1 = datetime.now()\n",
    "grpm_dataset['pmidmesh'] = grpm_dataset['pmids']+grpm_dataset['mesh']\n",
    "print(datetime.now() -time_1)\n",
    "grpm_dataset.head()\n",
    "\n",
    "time_1 = datetime.now()\n",
    "grpm_gene_metrics = grpm_dataset[['gene', 'rsid', 'pmids', 'mesh', 'pmidmesh']].groupby('gene').describe(include='all')\n",
    "print(datetime.now() -time_1)\n",
    "grpm_gene_metrics.head()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T14:36:00.860179700Z",
     "start_time": "2024-04-13T14:34:48.118711300Z"
    }
   },
   "execution_count": 23
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Filter GRPM dataset by MeSH [required]"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: total: 0 ns\n",
      "Wall time: 0 ns\n",
      "0:00:01.934377\n"
     ]
    }
   ],
   "source": [
    "%time\n",
    "\n",
    "# Preprocessing: \n",
    "time_1 = datetime.now()\n",
    "genes = grpm_dataset.gene.drop_duplicates().to_list()\n",
    "\n",
    "# Subset GRPM ds with ref-mesh\n",
    "mask_full = grpm_dataset['mesh'].isin(ref_mesh_list)\n",
    "grpm_match_full = grpm_dataset[mask_full].drop_duplicates()\n",
    "\n",
    "# save GRPM_gene\n",
    "grpm_match_full[['gene', 'rsid', 'pmids', 'mesh']].to_csv(directory+'/grpmx_filtered_output.csv')\n",
    "print(datetime.now() -time_1)\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T13:59:59.237754400Z",
     "start_time": "2024-04-13T13:59:57.299596100Z"
    }
   },
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Get GRPM Subset Metrics [required]"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "            gene         rsid     pmids                       mesh  \\\n45        MT-ND1  rs111033358  33468709  Diabetes Mellitus, Type 2   \n51        MT-ND1  rs111033358  33468709                    Insulin   \n52        MT-ND1  rs111033358  33468709         Insulin Resistance   \n1262      MT-ND1    rs1599988  28693754                Weight Gain   \n1281      MT-ND1   rs28358585  28693754                Weight Gain   \n...          ...          ...       ...                        ...   \n16608733     EVC    rs2291157  23139751                     Health   \n16608800     EVC    rs6414624  29273463  Diabetes Mellitus, Type 2   \n16609420    TMC3  rs150843673  29691411  Diabetes Mellitus, Type 2   \n16609428    TMC3  rs150843673  29691411                    Obesity   \n16609445    TMC3  rs150843673  33230300               Hypertension   \n\n                                   pmidmesh  \n45        33468709Diabetes Mellitus, Type 2  \n51                          33468709Insulin  \n52               33468709Insulin Resistance  \n1262                    28693754Weight Gain  \n1281                    28693754Weight Gain  \n...                                     ...  \n16608733                     23139751Health  \n16608800  29273463Diabetes Mellitus, Type 2  \n16609420  29691411Diabetes Mellitus, Type 2  \n16609428                    29691411Obesity  \n16609445               33230300Hypertension  \n\n[229421 rows x 5 columns]",
      "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>gene</th>\n      <th>rsid</th>\n      <th>pmids</th>\n      <th>mesh</th>\n      <th>pmidmesh</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>45</th>\n      <td>MT-ND1</td>\n      <td>rs111033358</td>\n      <td>33468709</td>\n      <td>Diabetes Mellitus, Type 2</td>\n      <td>33468709Diabetes Mellitus, Type 2</td>\n    </tr>\n    <tr>\n      <th>51</th>\n      <td>MT-ND1</td>\n      <td>rs111033358</td>\n      <td>33468709</td>\n      <td>Insulin</td>\n      <td>33468709Insulin</td>\n    </tr>\n    <tr>\n      <th>52</th>\n      <td>MT-ND1</td>\n      <td>rs111033358</td>\n      <td>33468709</td>\n      <td>Insulin Resistance</td>\n      <td>33468709Insulin Resistance</td>\n    </tr>\n    <tr>\n      <th>1262</th>\n      <td>MT-ND1</td>\n      <td>rs1599988</td>\n      <td>28693754</td>\n      <td>Weight Gain</td>\n      <td>28693754Weight Gain</td>\n    </tr>\n    <tr>\n      <th>1281</th>\n      <td>MT-ND1</td>\n      <td>rs28358585</td>\n      <td>28693754</td>\n      <td>Weight Gain</td>\n      <td>28693754Weight Gain</td>\n    </tr>\n    <tr>\n      <th>...</th>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n      <td>...</td>\n    </tr>\n    <tr>\n      <th>16608733</th>\n      <td>EVC</td>\n      <td>rs2291157</td>\n      <td>23139751</td>\n      <td>Health</td>\n      <td>23139751Health</td>\n    </tr>\n    <tr>\n      <th>16608800</th>\n      <td>EVC</td>\n      <td>rs6414624</td>\n      <td>29273463</td>\n      <td>Diabetes Mellitus, Type 2</td>\n      <td>29273463Diabetes Mellitus, Type 2</td>\n    </tr>\n    <tr>\n      <th>16609420</th>\n      <td>TMC3</td>\n      <td>rs150843673</td>\n      <td>29691411</td>\n      <td>Diabetes Mellitus, Type 2</td>\n      <td>29691411Diabetes Mellitus, Type 2</td>\n    </tr>\n    <tr>\n      <th>16609428</th>\n      <td>TMC3</td>\n      <td>rs150843673</td>\n      <td>29691411</td>\n      <td>Obesity</td>\n      <td>29691411Obesity</td>\n    </tr>\n    <tr>\n      <th>16609445</th>\n      <td>TMC3</td>\n      <td>rs150843673</td>\n      <td>33230300</td>\n      <td>Hypertension</td>\n      <td>33230300Hypertension</td>\n    </tr>\n  </tbody>\n</table>\n<p>229421 rows × 5 columns</p>\n</div>"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time\n",
    "grpm_match_full['pmidmesh'] = grpm_match_full['pmids']+grpm_match_full['mesh']\n",
    "time_1 = datetime.now()\n",
    "grpm_match_gene_metrics = grpm_match_full[['gene', 'rsid', 'pmids', 'mesh', 'pmidmesh']].groupby('gene').describe(include='all')\n",
    "print(datetime.now() -time_1)\n",
    "grpm_match_gene_metrics.head()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T14:00:13.072916600Z",
     "start_time": "2024-04-13T14:00:13.017576600Z"
    }
   },
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Other GRPM Subset Metrics [optional]"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:01:00.944193\n"
     ]
    }
   ],
   "source": [
    "time_1 = datetime.now()\n",
    "grpm_match_rsid_metrics = grpm_match_full[['rsid', 'pmids']].groupby('rsid').describe(include='all')#agg(lambda x: x.unique())\n",
    "print(datetime.now() -time_1)\n",
    "#grpm_match_rsid_metrics.head()\n",
    "grpm_match_rsid_metrics_sort = grpm_match_rsid_metrics.sort_values(by=('pmids','unique'), ascending= False)\n",
    "top10rsid = list(grpm_match_rsid_metrics_sort[:10].index)\n",
    "top10rsid"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T14:19:26.247705800Z",
     "start_time": "2024-04-13T14:18:25.293240400Z"
    }
   },
   "execution_count": 11
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:00:00.589811\n"
     ]
    },
    {
     "data": {
      "text/plain": "['Diabetes Mellitus, Type 2',\n 'Obesity',\n 'Inflammation',\n 'Body Mass Index',\n 'Hypertension',\n 'Insulin',\n 'Insulin Resistance',\n 'Family Health',\n 'Diet',\n 'Metabolic Syndrome']"
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "time_1 = datetime.now()\n",
    "grpm_match_mesh_metrics = grpm_match_full[['pmids', 'mesh']].groupby('mesh').describe(include='all')#agg(lambda x: x.unique())\n",
    "print(datetime.now() -time_1)\n",
    "#grpm_match_mesh_metrics.head()\n",
    "\n",
    "grpm_match_mesh_metrics_sort = grpm_match_mesh_metrics.sort_values(by=('pmids','unique'), ascending= False)\n",
    "top10mesh = list(grpm_match_mesh_metrics_sort[:10].index)\n",
    "top10mesh"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T14:27:36.864074500Z",
     "start_time": "2024-04-13T14:27:36.279787100Z"
    }
   },
   "execution_count": 21
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "time_1 = datetime.now()\n",
    "grpm_match_pmids_metrics = grpm_match_full[['gene', 'rsid', 'pmids', 'mesh']].groupby('pmids').describe(include='all')#agg(lambda x: x.unique())\n",
    "print(datetime.now() -time_1)\n",
    "grpm_match_pmids_metrics.head()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "0:00:57.947976 x 3\n",
    "0:01:14.471060 x 4"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Run Survey"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 15519/15519 [10:28<00:00, 24.67it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:00:00.022980\n",
      "time batch: 0:10:30.102083\n",
      "runtime/gene: 0:00:00.040602\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "#---------------------------------------------\n",
    "# Edit saving options:\n",
    "checkpoint = 500 # write Report each \"n\" genes\n",
    "\n",
    "include_top10 = True # for a faster job\n",
    "save_plot = False  # only if include_top10 = True\n",
    "\n",
    "run_sample = False # set True just to run a test\n",
    "sample_size = 10\n",
    "\n",
    "df_report_complete = pd.DataFrame()\n",
    "#---------------------------------------------\n",
    "\n",
    "\n",
    "# Define grpm subset\n",
    "if 'subset_grpm' in locals() and subset_grpm:\n",
    "    GRPM_report = GRPM_report_subset\n",
    "    grpm_dataset = grpm_subset\n",
    "\n",
    "# define gene list\n",
    "import random\n",
    "if run_sample:\n",
    "    genes = random.sample(grpm_genes_list[:], sample_size)\n",
    "else:\n",
    "    genes = grpm_genes_list[gene_start:len(grpm_genes_list)]\n",
    "\n",
    "time_start = datetime.now()\n",
    "\n",
    "for gene in tqdm(genes):\n",
    "\n",
    "    time_alpha = datetime.now()\n",
    "    timestamp = time_alpha.strftime('%Y%m%d%H%M%S')\n",
    "\n",
    "    # Pre-Selection Metrics ===========\n",
    "    grpm_gene_metrics_gene = grpm_gene_metrics.loc[gene]\n",
    "    starting_pmid          =  grpm_gene_metrics_gene.pmids.loc['unique']\n",
    "    starting_mesh          =  grpm_gene_metrics_gene.mesh.loc['unique']\n",
    "    lit1_rsid              =  grpm_gene_metrics_gene.rsid.loc['unique']\n",
    "    starting_pmidmesh      =  grpm_gene_metrics_gene.pmidmesh.loc['unique']\n",
    "    \n",
    "    #  Post-Selection Metrics ===========\n",
    "    if gene in grpm_match_gene_metrics.index:\n",
    "        grpm_match_gene_metrics_gene = grpm_match_gene_metrics.loc[gene]\n",
    "        matching_rsid                = grpm_match_gene_metrics_gene.rsid.loc['unique']\n",
    "        dropped_rsid                 = lit1_rsid - matching_rsid\n",
    "        matching_pmids               = grpm_match_gene_metrics_gene.pmids.loc['unique']\n",
    "        matching_mesh                = grpm_match_gene_metrics_gene.mesh.loc['unique']\n",
    "        matching_pmidmesh            = grpm_match_gene_metrics_gene.pmidmesh.loc['unique']\n",
    "    else:\n",
    "        matching_rsid                = 0\n",
    "        dropped_rsid                 = lit1_rsid - 0\n",
    "        matching_pmids               = 0\n",
    "        matching_mesh                = 0\n",
    "        matching_pmidmesh            = 0\n",
    "                \n",
    "\n",
    "    #------------------------\n",
    "    if not include_top10:\n",
    "        matching_rsid_pmid10  = 'na'\n",
    "        matching_rsid_pmid100 = 'na'\n",
    "        top10rsid             = 'na'\n",
    "        top10mesh             = 'na'\n",
    "    else:\n",
    "        #=====================================\n",
    "        # Get gene_grpm (slow step)\n",
    "        grpm_gene = grpm_match_full.loc[grpm_match_full['gene'] == gene]\n",
    "        #print(datetime.now() -time_alpha)\n",
    "    \n",
    "        dfmatch_full = grpm_gene\n",
    "        dfmatch_less = dfmatch_full[['pmids', 'rsid', 'mesh']].drop_duplicates()\n",
    "        #=====================================\n",
    "\n",
    "        #=====================================\n",
    "        ## 1. groupby.describe analysis by [rsid]\n",
    "        dfmatch_less_rsid = dfmatch_less.groupby('rsid').describe().reset_index()\n",
    "        dfmatch_less_rsid.columns = dfmatch_less_rsid.columns.to_flat_index()\n",
    "        new_column_names = ['rsid', 'pmid-count', 'pmid-unique','pmid-top','pmid-freq','mesh-count', 'mesh-unique','mesh-top','mesh-freq']\n",
    "        dfmatch_less_rsid.columns = new_column_names\n",
    "\n",
    "        ### statistics:\n",
    "        #grpm_match_rsid_metrics_gene = grpm_match_rsid_metrics.loc[gene]\n",
    "        matching_rsid_pmid10 = len(dfmatch_less_rsid[dfmatch_less_rsid['pmid-unique']>10])\n",
    "        matching_rsid_pmid100 = len(dfmatch_less_rsid[dfmatch_less_rsid['pmid-unique']>100])\n",
    "\n",
    "        ### sorting, top10\n",
    "        dfmatch_less_rsidless = dfmatch_less_rsid[['rsid','pmid-unique','mesh-unique']]\n",
    "        dfmatch_less_rsidlesssort = dfmatch_less_rsidless.sort_values(by='pmid-unique', ascending= False).reset_index(drop=True)\n",
    "        top10rsid = dfmatch_less_rsidlesssort['rsid'][:10].tolist()\n",
    "        #==================================\n",
    "        \n",
    "        #==================================\n",
    "        ## 2. groupby.describe analysis by [mesh]\n",
    "        dfmatch_less_mesh = dfmatch_less.groupby('mesh').describe().reset_index()\n",
    "        dfmatch_less_mesh.columns = dfmatch_less_mesh.columns.to_flat_index()\n",
    "        # simplify df.groupby.describe, convert Multicolumn to single column \n",
    "        new_column_names = ['mesh', 'pmid-count', 'pmid-unique','pmid-top','pmid-freq','rsid-count', 'rsid-unique','rsid-top','rsid-freq']\n",
    "        dfmatch_less_mesh.columns = new_column_names\n",
    "        dfmatch_less_mesh_less = dfmatch_less_mesh[['mesh','pmid-unique','rsid-unique']]\n",
    "        \n",
    "        ### add frequency, top10\n",
    "        samplepmid_count = len(dfmatch_less.pmids.drop_duplicates())\n",
    "        dfmatch_less_mesh_less_frq = dfmatch_less_mesh_less.copy()\n",
    "        mesh_frq = dfmatch_less_mesh_less_frq.loc[:,'pmid-unique'].astype(float)/samplepmid_count\n",
    "        dfmatch_less_mesh_less_frq.loc[:,'mesh frequency'] = round(mesh_frq,3)#*100\n",
    "        dfmatch_less_mesh_less_frqsort = dfmatch_less_mesh_less_frq.sort_values(by='pmid-unique',ascending=False).reset_index(drop=True)\n",
    "        top10mesh = dfmatch_less_mesh_less_frqsort['mesh'][:10].tolist()\n",
    "        #==================================\n",
    "\n",
    "        if save_plot:\n",
    "            # create a scatter plot\n",
    "            x = dfmatch_less_mesh_less_frqsort['mesh'].head(30)\n",
    "            y = dfmatch_less_mesh_less_frqsort['pmid-unique'].head(30)\n",
    "            plt.figure(figsize=(5, 8))\n",
    "            plt.title('Scatter Plot: '+gene+' pmid-mesh (filtered)', loc='center',pad=10)\n",
    "            plt.scatter(y, x)\n",
    "            plt.gca().invert_yaxis()\n",
    "            #plt.subplots_adjust(left=0.3, right=0.9, bottom=0.3, top=0.9)\n",
    "            #plt.xticks(rotation=90)\n",
    "            plt.tick_params(axis='x', which='both', top=True, bottom=False, labeltop=True, labelbottom=False)\n",
    "            plt.xlabel('pmid count', position=(0.5, 1.08))\n",
    "            ax = plt.gca()\n",
    "            ax.xaxis.set_label_position('top')\n",
    "            #plt.show()\n",
    "            plt.savefig(directory+'/'+gene+'_mesh_plot_'+timestamp+'_filtered.png',dpi=120, bbox_inches = \"tight\")\n",
    "            plt.close()\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "\n",
    "    # Collect REPORT data ==================================\n",
    "\n",
    "    report = { 'reference_mesh': ref_mesh_n,\n",
    "               'starting_pmidmesh': starting_pmidmesh,\n",
    "               'starting_pmid' : starting_pmid,\n",
    "               'starting_mesh': starting_mesh,\n",
    "               'starting_rsid': lit1_rsid,\n",
    "               'matching_pmidmesh': matching_pmidmesh,\n",
    "               'matching_pmids': matching_pmids,\n",
    "               'matching_mesh': matching_mesh,\n",
    "               'matching_rsid': matching_rsid,\n",
    "               'dropped_rsid': dropped_rsid,\n",
    "               'matching_mesh_ratio': round((matching_mesh/starting_mesh),3),\n",
    "               'matching_pmids_ratio': round((matching_pmids/starting_pmid),3),\n",
    "               'matching_pmidmesh_ratio':  round((matching_pmidmesh/starting_pmidmesh),3),\n",
    "               'matching_rsid_ratio': round((matching_rsid/lit1_rsid),3),\n",
    "               'matching_rsid_pmid10': matching_rsid_pmid10,\n",
    "               'matching_rsid_pmid100': matching_rsid_pmid100,\n",
    "               'matching_top10mesh':str(top10mesh),\n",
    "               'matching_top10rsid':str(top10rsid),\n",
    "               }\n",
    "    \n",
    "    # UPDATE REPORT ------------\n",
    "    df_new_raw = pd.DataFrame(report, index=[gene])\n",
    "    df_report_complete = pd.concat([df_report_complete, df_new_raw])\n",
    "\n",
    "    full_runtime = datetime.now() - time_alpha\n",
    "    #print((gene+'_runtime:').ljust(18)+ str(full_runtime).ljust(15), ' Genes processed:', genes.index(gene), 'on', len(genes))\n",
    "    total_seconds = full_runtime.total_seconds()\n",
    "\n",
    "    # save checkpoint----------------------\n",
    "    if genes.index(gene) > 1 and genes.index(gene) % checkpoint == 0:\n",
    "        df_report_complete.to_csv(directory+'/GRPMX_report.csv')\n",
    "        #print(\"saved checkpoint\")\n",
    "    else:\n",
    "        pass\n",
    "    #==================================\n",
    "\n",
    "\n",
    "# Save report csv  (saving translate version is a code atavism)\n",
    "df_report_complete.T.to_csv(directory+'/GRPMX_report.csv')\n",
    "\n",
    "# #Update gene values (remove previous gene entry)\n",
    "GRPMX_report = pd.read_csv(directory+'/GRPMX_report.csv', index_col=0)\n",
    "time_load_1 = datetime.now()\n",
    "for gene in grpm_genes_list:\n",
    "    if gene+'.1' in GRPMX_report.columns:\n",
    "        GRPMX_report = GRPMX_report.drop(columns = gene)\n",
    "        GRPMX_report = GRPMX_report.rename(columns={gene+'.1': gene})\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "print(datetime.now() - time_load_1)\n",
    "GRPMX_report.to_csv(directory+'/GRPMX_report.csv')\n",
    "\n",
    "time_finish = datetime.now()\n",
    "time_batch = time_finish - time_start\n",
    "\n",
    "if os.path.isfile('run_time.txt'):\n",
    "    with open('run_time.txt', 'a') as file:\n",
    "        file.write(topic_tag+':\\n\\ttime batch: '+str(time_batch)+'\\n\\truntime/gene: '+str(time_batch/len(genes))+'\\n\\n')\n",
    "else:\n",
    "    with open('run_time.txt', 'w') as file:\n",
    "        file.write(topic_tag+':\\n\\ttime batch: '+str(time_batch)+'\\n\\truntime/gene: '+str(time_batch/len(genes))+'\\n\\n')\n",
    "\n",
    "print('time batch:',time_batch)\n",
    "print('runtime/gene:', time_batch/len(genes))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T15:38:37.459805700Z",
     "start_time": "2024-04-13T15:28:07.338376900Z"
    }
   },
   "execution_count": 61
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "         reference_mesh  starting_pmidmesh  starting_pmid  starting_mesh  \\\nMT-ND1              243               3603            269            764   \nMT-ND2              243               3958            292            795   \nMT-CO1              243               2308            163            558   \nMT-CO2              243               1330             94            422   \nMT-ATP8             243               1802            127            518   \n...                 ...                ...            ...            ...   \nCRISP1              243                 22              2             19   \nTMC3                243                 50              4             40   \nMPIG6B              243                210             16            107   \nPDP1                243                241             19            179   \nTAS2R60             243                 85              6             56   \n\n         starting_rsid  matching_pmidmesh  matching_pmids  matching_mesh  \\\nMT-ND1             105                 51              40             40   \nMT-ND2             123                 57              46             46   \nMT-CO1              97                 29              26             26   \nMT-CO2              67                 21              17             17   \nMT-ATP8             77                 23              19             19   \n...                ...                ...             ...            ...   \nCRISP1               2                  0               0              0   \nTMC3                 3                  3               2              2   \nMPIG6B               7                  3               3              3   \nPDP1                13                  2               2              2   \nTAS2R60              3                  4               2              2   \n\n         matching_rsid  dropped_rsid  matching_mesh_ratio  \\\nMT-ND1              27            78                0.052   \nMT-ND2              29            94                0.058   \nMT-CO1              22            75                0.047   \nMT-CO2              14            53                0.040   \nMT-ATP8             14            63                0.037   \n...                ...           ...                  ...   \nCRISP1               0             0                0.000   \nTMC3                 1             2                0.050   \nMPIG6B               4             3                0.028   \nPDP1                 1            12                0.011   \nTAS2R60              1             2                0.036   \n\n         matching_pmids_ratio  matching_pmidmesh_ratio  matching_rsid_ratio  \\\nMT-ND1                  0.149                    0.014                0.257   \nMT-ND2                  0.158                    0.014                0.236   \nMT-CO1                  0.160                    0.013                0.227   \nMT-CO2                  0.181                    0.016                0.209   \nMT-ATP8                 0.150                    0.013                0.182   \n...                       ...                      ...                  ...   \nCRISP1                  0.000                    0.000                0.000   \nTMC3                    0.500                    0.060                0.333   \nMPIG6B                  0.188                    0.014                0.571   \nPDP1                    0.105                    0.008                0.077   \nTAS2R60                 0.333                    0.047                0.333   \n\n        matching_rsid_pmid10 matching_rsid_pmid100 matching_top10mesh  \\\nMT-ND1                    na                    na                 na   \nMT-ND2                    na                    na                 na   \nMT-CO1                    na                    na                 na   \nMT-CO2                    na                    na                 na   \nMT-ATP8                   na                    na                 na   \n...                      ...                   ...                ...   \nCRISP1                    na                    na                 na   \nTMC3                      na                    na                 na   \nMPIG6B                    na                    na                 na   \nPDP1                      na                    na                 na   \nTAS2R60                   na                    na                 na   \n\n        matching_top10rsid  \nMT-ND1                  na  \nMT-ND2                  na  \nMT-CO1                  na  \nMT-CO2                  na  \nMT-ATP8                 na  \n...                    ...  \nCRISP1                  na  \nTMC3                    na  \nMPIG6B                  na  \nPDP1                    na  \nTAS2R60                 na  \n\n[15519 rows x 18 columns]",
      "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>reference_mesh</th>\n      <th>starting_pmidmesh</th>\n      <th>starting_pmid</th>\n      <th>starting_mesh</th>\n      <th>starting_rsid</th>\n      <th>matching_pmidmesh</th>\n      <th>matching_pmids</th>\n      <th>matching_mesh</th>\n      <th>matching_rsid</th>\n      <th>dropped_rsid</th>\n      <th>matching_mesh_ratio</th>\n      <th>matching_pmids_ratio</th>\n      <th>matching_pmidmesh_ratio</th>\n      <th>matching_rsid_ratio</th>\n      <th>matching_rsid_pmid10</th>\n      <th>matching_rsid_pmid100</th>\n      <th>matching_top10mesh</th>\n      <th>matching_top10rsid</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>MT-ND1</th>\n      <td>243</td>\n      <td>3603</td>\n      <td>269</td>\n      <td>764</td>\n      <td>105</td>\n      <td>51</td>\n      <td>40</td>\n      <td>40</td>\n      <td>27</td>\n      <td>78</td>\n      <td>0.052</td>\n      <td>0.149</td>\n      <td>0.014</td>\n      <td>0.257</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>MT-ND2</th>\n      <td>243</td>\n      <td>3958</td>\n      <td>292</td>\n      <td>795</td>\n      <td>123</td>\n      <td>57</td>\n      <td>46</td>\n      <td>46</td>\n      <td>29</td>\n      <td>94</td>\n      <td>0.058</td>\n      <td>0.158</td>\n      <td>0.014</td>\n      <td>0.236</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>MT-CO1</th>\n      <td>243</td>\n      <td>2308</td>\n      <td>163</td>\n      <td>558</td>\n      <td>97</td>\n      <td>29</td>\n      <td>26</td>\n      <td>26</td>\n      <td>22</td>\n      <td>75</td>\n      <td>0.047</td>\n      <td>0.160</td>\n      <td>0.013</td>\n      <td>0.227</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>MT-CO2</th>\n      <td>243</td>\n      <td>1330</td>\n      <td>94</td>\n      <td>422</td>\n      <td>67</td>\n      <td>21</td>\n      <td>17</td>\n      <td>17</td>\n      <td>14</td>\n      <td>53</td>\n      <td>0.040</td>\n      <td>0.181</td>\n      <td>0.016</td>\n      <td>0.209</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>MT-ATP8</th>\n      <td>243</td>\n      <td>1802</td>\n      <td>127</td>\n      <td>518</td>\n      <td>77</td>\n      <td>23</td>\n      <td>19</td>\n      <td>19</td>\n      <td>14</td>\n      <td>63</td>\n      <td>0.037</td>\n      <td>0.150</td>\n      <td>0.013</td>\n      <td>0.182</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</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    </tr>\n    <tr>\n      <th>CRISP1</th>\n      <td>243</td>\n      <td>22</td>\n      <td>2</td>\n      <td>19</td>\n      <td>2</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0</td>\n      <td>0.000</td>\n      <td>0.000</td>\n      <td>0.000</td>\n      <td>0.000</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>TMC3</th>\n      <td>243</td>\n      <td>50</td>\n      <td>4</td>\n      <td>40</td>\n      <td>3</td>\n      <td>3</td>\n      <td>2</td>\n      <td>2</td>\n      <td>1</td>\n      <td>2</td>\n      <td>0.050</td>\n      <td>0.500</td>\n      <td>0.060</td>\n      <td>0.333</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>MPIG6B</th>\n      <td>243</td>\n      <td>210</td>\n      <td>16</td>\n      <td>107</td>\n      <td>7</td>\n      <td>3</td>\n      <td>3</td>\n      <td>3</td>\n      <td>4</td>\n      <td>3</td>\n      <td>0.028</td>\n      <td>0.188</td>\n      <td>0.014</td>\n      <td>0.571</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>PDP1</th>\n      <td>243</td>\n      <td>241</td>\n      <td>19</td>\n      <td>179</td>\n      <td>13</td>\n      <td>2</td>\n      <td>2</td>\n      <td>2</td>\n      <td>1</td>\n      <td>12</td>\n      <td>0.011</td>\n      <td>0.105</td>\n      <td>0.008</td>\n      <td>0.077</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n    <tr>\n      <th>TAS2R60</th>\n      <td>243</td>\n      <td>85</td>\n      <td>6</td>\n      <td>56</td>\n      <td>3</td>\n      <td>4</td>\n      <td>2</td>\n      <td>2</td>\n      <td>1</td>\n      <td>2</td>\n      <td>0.036</td>\n      <td>0.333</td>\n      <td>0.047</td>\n      <td>0.333</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n      <td>na</td>\n    </tr>\n  </tbody>\n</table>\n<p>15519 rows × 18 columns</p>\n</div>"
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_report_complete"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-04-13T15:02:36.707000600Z",
     "start_time": "2024-04-13T15:02:36.707000600Z"
    }
   },
   "execution_count": 56
  },
  {
   "cell_type": "markdown",
   "source": [
    "GRPM-Survey Version 2.0\n",
    "\n",
    "OB-BMI (full)/\n",
    "time batch: 0:42:36.831920\n",
    "runtime/gene: 0:00:00.164755\n",
    "\n",
    "======================================\n",
    "GRPM-Survey Version 3.0\n",
    "OB-BMI (partial)\n",
    "time batch: 0:00:49.547031\n",
    "runtime/gene: 0:00:00.003193\n",
    "\n",
    "OB-BMI (full)\n",
    "time batch: 0:10:30.102083\n",
    "runtime/gene: 0:00:00.040602\\78/"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "pd.read_csv(directory+'/GRPMX_report.csv').T.head(20)"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_read = pd.read_csv(directory+'/grpmx_filtered_output.csv', index_col=0)\n",
    "print('genes matching:', df_read.gene.nunique())\n",
    "print('mesh matching:', df_read.mesh.nunique())\n",
    "print('apply threshold in Analyzer Module')\n",
    "df_read"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Visualize GRPMX_report.csv\n",
    "GRPMX_report = pd.read_csv(directory+'/GRPMX_report.csv', index_col=0).transpose().reset_index().rename(columns={'index':'gene'})\n",
    "GRPMX_report.gene.drop_duplicates().to_clipboard()\n",
    "print('Genes matching:',len(GRPMX_report.gene.drop_duplicates()))\n",
    "\n",
    "GRPMX_report[['reference_mesh', 'starting_pmidmesh', 'starting_pmid','starting_mesh','starting_rsid', 'matching_pmidmesh', 'matching_pmids', 'matching_mesh','matching_rsid', 'dropped_rsid']] = GRPMX_report[['reference_mesh', 'starting_pmidmesh', 'starting_pmid','starting_mesh','starting_rsid', 'matching_pmidmesh', 'matching_pmids', 'matching_mesh','matching_rsid', 'dropped_rsid']].astype(int)\n",
    "\n",
    "GRPMX_report[['matching_mesh_ratio', 'matching_pmids_ratio','matching_pmidmesh_ratio', 'matching_rsid_ratio']] = GRPMX_report[['matching_mesh_ratio', 'matching_pmids_ratio','matching_pmidmesh_ratio','matching_rsid_ratio']].astype(float)\n",
    "\n",
    "columns_to_keep = ['matching_pmids','matching_pmids_ratio','matching_mesh','matching_rsid']\n",
    "GRPMX_report_less = GRPMX_report[columns_to_keep]\n",
    "\n",
    "sorting_column = 'matching_pmids'\n",
    "GRPMX_report_sort = GRPMX_report.sort_values(by=sorting_column, ascending=False)\n",
    "\n",
    "columns_to_display = ['gene', 'matching_pmidmesh', 'matching_pmids',\n",
    "                      'matching_mesh', 'matching_rsid', 'dropped_rsid', 'matching_mesh_ratio',\n",
    "                      'matching_pmids_ratio', 'matching_pmidmesh_ratio',\n",
    "                      'matching_rsid_ratio']\n",
    "GRPMX_report_display = GRPMX_report[columns_to_display]\n",
    "GRPMX_report_display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Matching PMIDs in Database\n",
    "GRPMX_report_sort = GRPMX_report.sort_values(by= 'matching_pmids',ascending=False)\n",
    "\n",
    "x = GRPMX_report_sort.gene.iloc[:40]\n",
    "y = GRPMX_report_sort['matching_pmids'].iloc[:40]\n",
    "plt.figure(figsize=(5, len(x)*0.2))\n",
    "plt.title('Matching PMIDs in Database', loc='center',pad=10)\n",
    "\n",
    "plt.barh(x,y)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.tick_params(axis='x', which='both', top=True, bottom=False, labeltop=True, labelbottom=False)\n",
    "#plt.xlabel('pmid count', position=(0.5, 1.08))\n",
    "plt.ylabel('genes')\n",
    "plt.xlabel('matching pmid', position=(0.5, 1.08))\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_label_position('top')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Add \"interest value\" to report:----------------------------------------------------------\n",
    "max_match_pmids = int(GRPMX_report['matching_pmids'].max())\n",
    "GRPMX_report_int = GRPMX_report\n",
    "GRPMX_report_int['matching_pmids_score'] = round((GRPMX_report_int['matching_pmids']/max_match_pmids),3)\n",
    "\n",
    "GRPMX_report_int['interest_value'] = round(GRPMX_report_int['matching_pmids_score'] * GRPMX_report_int['matching_pmids_ratio'],3)\n",
    "\n",
    "GRPMX_report_int.set_index('gene').sort_values(by='interest_value')#.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Matching PMIDs in Database\n",
    "GRPMX_report_sort = GRPMX_report.sort_values(by= 'matching_pmids_index',ascending=False)\n",
    "\n",
    "x = GRPMX_report_sort.gene.iloc[:100]\n",
    "y = GRPMX_report_sort['matching_pmids_index'].iloc[:100]\n",
    "plt.figure(figsize=(5, len(x)*0.2))\n",
    "plt.title('Matching PMIDs in Database', loc='center',pad=10)\n",
    "\n",
    "plt.barh(x,y)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.tick_params(axis='x', which='both', top=True, bottom=False, labeltop=True, labelbottom=False)\n",
    "#plt.xlabel('pmid count', position=(0.5, 1.08))\n",
    "plt.ylabel('genes')\n",
    "plt.xlabel('matching pmid', position=(0.5, 1.08))\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_label_position('top')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Extra"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Simple GRPM Subsetting"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "grpm_dataset.head()"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## filter by mesh"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "# filtering source dataset\n",
    "import time\n",
    "timea = time.time()\n",
    "\n",
    "my_mesh = 'Heart Failure'\n",
    "#nbib_subset = pd.DataFrame(columns= nbib_dataset.columns)\n",
    "filteres_grpm = grpm_dataset[grpm_dataset.mesh == my_mesh].reset_index(drop=True)\n",
    "\n",
    "print((time.time()-timea)/60,'minutes')\n",
    "filteres_grpm.to_csv('filteres_grpm_heart_fail.csv') #= pd.read_csv(filteres_grpm)\n",
    "filteres_grpm"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "# import LitVat-PubMed Dataset (GRPM)\n",
    "filteres_grpm = pd.read_csv('filteres_grpm_heart_fail.csv', index_col=0)\n",
    "filteres_grpm.pmids = filteres_grpm.pmids.astype('str')  # convert PMIDs to str\n",
    "filteres_grpm"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "#Analyze data with \"groupby.describe\" method\n",
    "\n",
    "## 1. groupby.describe analysis by [pmids]\n",
    "filteres_grpm_gene = filteres_grpm.groupby('gene').describe().reset_index()#.reset_index(drop=True)\n",
    "filteres_grpm_gene[['gene','pmids']].sort_values(by=('pmids', 'unique'), ascending=False).reset_index(drop=True)"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## filter by gene"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "# filtering source dataset\n",
    "import time\n",
    "timea = time.time()\n",
    "my_mesh = ref.mesh[0]\n",
    "my_gene = 'GLA'\n",
    "#nbib_subset = pd.DataFrame(columns= nbib_dataset.columns)\n",
    "filteres_grpm = grpm_dataset[grpm_dataset.gene == my_gene].reset_index(drop=True)\n",
    "print((time.time()-timea)/60,'minutes')\n",
    "\n",
    "filteres_grpm.to_csv('filteres_grpm.csv') #= pd.read_csv(filteres_grpm)"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "# import LitVat-PubMed Dataset (GRPM)\n",
    "filteres_grpm = pd.read_csv('filteres_grpm.csv', index_col=0)\n",
    "filteres_grpm_sub = filteres_grpm[filteres_grpm.mesh.str.contains('Fabry')].reset_index(drop=True)\n",
    "filteres_grpm_sub.pmids = filteres_grpm.pmids.astype('str')  # convert PMIDs to str\n",
    "filteres_grpm_sub"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "#Analyze data with \"groupby.describe\" method\n",
    "\n",
    "## 1. groupby.describe analysis by [pmids]\n",
    "filteres_grpm_sub_pmids = filteres_grpm.groupby('pmids').describe().reset_index()\n",
    "filteres_grpm_sub_pmids"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "## 1. groupby.describe analysis by [rsid]\n",
    "filteres_grpm_sub_rsid = filteres_grpm.groupby('rsid').describe().reset_index()\n",
    "filteres_grpm_sub_rsid"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "## 1. groupby.describe analysis by [rsid]\n",
    "filteres_grpm_mesh = filteres_grpm.groupby('mesh').describe().reset_index()\n",
    "filteres_grpm_mesh"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "filteres_grpm_[filteres_grpm_.mesh.str.contains('Fabry')]"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "\n",
    "filteres_grpm_.columns = filteres_grpm_.columns.to_flat_index()\n",
    "#new_column_names = ['rsid', 'pmid-count', 'pmid-unique','pmid-top','pmid-freq','mesh-count', 'mesh-unique','mesh-top','mesh-freq']\n",
    "filteres_grpm_.columns = filteres_grpm_\n",
    "#------------------\n",
    "filteres_grpm_"
   ],
   "metadata": {
    "collapsed": false
   },
   "execution_count": null
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "import json\n",
    "pmid_list =[]\n",
    "\n",
    "for  pmid_id  in myjson.elemento[0]:\n",
    "    pmid_list.append(pmid_id)\n"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}