examples/local_hmmer_example.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example notebook tutorial\n",
    "## local HMMER walkthrough"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# system dependecies\n",
    "\n",
    "\n",
    "# library dependencies\n",
    "import duckdb as ddb\n",
    "import pandas as pd\n",
    "from joblib import Parallel, delayed\n",
    "from tqdm import tqdm\n",
    "\n",
    "# local dependencies/utils\n",
    "from pairpro.hmmer import local_hmmer_wrapper_example"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook walkthrough, we are assuming the user has downloaded pfam-A HMMs, which be downloaded by running the following command after clonning the repository:\n",
    "\n",
    "```bash\n",
    "python scripts/download_pfam.py\n",
    "```\n",
    "\n",
    "You will only need to run this command once, and you will find the HMMs in the `data` directory.\n",
    "\n",
    "We will also use an adapated version of fastest local HMMER wrapper, which assumes the input is a database of protein pairs, but we have other versions of the wrapper that can be used for other purposes. For more information, please refer to the Read the Docs page.\n",
    "\n",
    "The aim of this walkthrough is mainly to elucidate how the local HMMER wrapper works, and how to use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pfam path\n",
    "PFAM_PATH = \"../data/pfam/\"\n",
    "\n",
    "# any user-defined protein db\n",
    "SAMPLE_DB_PATH = \"../data/pairpro_50k.db\"\n",
    "\n",
    "# output path\n",
    "HMMER_OUTPUT_DIR = '../data/protein_pairs/'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data and preprocess it a bit"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "let's see how the sample database looks like. Let's us first make a connection to the database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn = ddb.connect(SAMPLE_DB_PATH, read_only=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make to a dataframe to see how the data looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_query = conn.execute(\"SELECT * FROM pairpro.proteins LIMIT 10\").df()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pid</th>\n",
       "      <th>taxid</th>\n",
       "      <th>pdb_id</th>\n",
       "      <th>alphafold_id</th>\n",
       "      <th>proteome</th>\n",
       "      <th>protein_seq</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>A0A6L5BYG5</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BYG5</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MLQRYLWKLLPKQQRAFLLGRLSVVDRQVVNKSMSANLQFPSSFAQ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>A0A6L5BTN2</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BTN2</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MEQQEAWQVLIVEDDQRLAELTRDYLEANGLRVAIEGNGALAAARI...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A0A6L5BN44</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BN44</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MQPFVIAPSILSADFARLGEEVDNVLAAGADFVHFDVMDNHYVPNL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>A0A6L5BWT7</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BWT7</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MQVESRPDKKSGRFFMRIGHGYDVHRFAEGDFITLGGVRIAHGFGL...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>A0A6L5BWU3</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BWU3</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MRPSEWFEGLRKIDINDLDTNNIGSWPPAIKALAGILLMVLVLGLG...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>A0A6L5BX00</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BX00</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MTPSLLMAVLASGFIYGITPGPGVLAVFGIGAARGRRAGAGFLCGH...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>A0A6L5C1M3</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5C1M3</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MSRLKNKYALITGGTSGIGLETARQFLAQGATVAITGRSESALAAA...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>A0A6L5BNJ3</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BNJ3</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MEQTKRVLVVEDDLHIADLICLHLRDEQFEVVHCADGDEGMRLLQQ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>A0A6L5BR68</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BR68</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MFTKQRLIIVATAVALLSGCASPNPYDNQGQADGGSQGMSKTAKYG...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>A0A6L5BVP2</td>\n",
       "      <td>104087</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A0A6L5BVP2</td>\n",
       "      <td>UP000475265</td>\n",
       "      <td>MQNPQNLIWIDLEMTGLNPDTDVIIEMATIVTDSDLNTLAEGPVIA...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          pid   taxid  pdb_id alphafold_id     proteome  \\\n",
       "0  A0A6L5BYG5  104087     NaN   A0A6L5BYG5  UP000475265   \n",
       "1  A0A6L5BTN2  104087     NaN   A0A6L5BTN2  UP000475265   \n",
       "2  A0A6L5BN44  104087     NaN   A0A6L5BN44  UP000475265   \n",
       "3  A0A6L5BWT7  104087     NaN   A0A6L5BWT7  UP000475265   \n",
       "4  A0A6L5BWU3  104087     NaN   A0A6L5BWU3  UP000475265   \n",
       "5  A0A6L5BX00  104087     NaN   A0A6L5BX00  UP000475265   \n",
       "6  A0A6L5C1M3  104087     NaN   A0A6L5C1M3  UP000475265   \n",
       "7  A0A6L5BNJ3  104087     NaN   A0A6L5BNJ3  UP000475265   \n",
       "8  A0A6L5BR68  104087     NaN   A0A6L5BR68  UP000475265   \n",
       "9  A0A6L5BVP2  104087     NaN   A0A6L5BVP2  UP000475265   \n",
       "\n",
       "                                         protein_seq  \n",
       "0  MLQRYLWKLLPKQQRAFLLGRLSVVDRQVVNKSMSANLQFPSSFAQ...  \n",
       "1  MEQQEAWQVLIVEDDQRLAELTRDYLEANGLRVAIEGNGALAAARI...  \n",
       "2  MQPFVIAPSILSADFARLGEEVDNVLAAGADFVHFDVMDNHYVPNL...  \n",
       "3  MQVESRPDKKSGRFFMRIGHGYDVHRFAEGDFITLGGVRIAHGFGL...  \n",
       "4  MRPSEWFEGLRKIDINDLDTNNIGSWPPAIKALAGILLMVLVLGLG...  \n",
       "5  MTPSLLMAVLASGFIYGITPGPGVLAVFGIGAARGRRAGAGFLCGH...  \n",
       "6  MSRLKNKYALITGGTSGIGLETARQFLAQGATVAITGRSESALAAA...  \n",
       "7  MEQTKRVLVVEDDLHIADLICLHLRDEQFEVVHCADGDEGMRLLQQ...  \n",
       "8  MFTKQRLIIVATAVALLSGCASPNPYDNQGQADGGSQGMSKTAKYG...  \n",
       "9  MQNPQNLIWIDLEMTGLNPDTDVIIEMATIVTDSDLNTLAEGPVIA...  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_query"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a pretty feature-rich dataframe, but we only read sequence pairs and their indexes for the hmmer_wrapper, so let's do a bit of preprocessing before we run HMMER against pfam!\n",
    "The wrapper expects a dataframe with the following columns:\n",
    "\n",
    "pid, protein_seq \n",
    "\n",
    "Where the pid is the protein ID, and protein_seq is the protein sequence. These two columns are only for proteins from pairs, so HMMER runs the most efficently, and we can get the most out of the wrapper. The wrapper will take a chunked list of PIDs, and query the databse for the sequences, and then run HMMER against the sequences via pfam."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins_in_pair_pids = conn.execute(f\"SELECT pid FROM pairpro.proteins LIMIT 4000\").df()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "chunk_size = 1000\n",
    "# chunking the PID so the worker function queries\n",
    "protein_pair_pid_chunks = [proteins_in_pair_pids[i:i + chunk_size]\n",
    "                        for i in range(0, len(proteins_in_pair_pids), chunk_size)]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "let's run the wrapper in parallel, and see how it works!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "njobs = 4 # number of jobs to run in parallel\n",
    "conn.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 25%|██▌       | 1/4 [01:08<03:25, 68.62s/it]\n"
     ]
    }
   ],
   "source": [
    "with tqdm(total=len(protein_pair_pid_chunks)) as pbar:\n",
    "        pbar.update(1)\n",
    "        Parallel(\n",
    "            n_jobs=njobs)(\n",
    "            delayed(local_hmmer_wrapper_example)(\n",
    "                chunk_index,\n",
    "                SAMPLE_DB_PATH,\n",
    "                protein_pair_pid_chunks,\n",
    "                PFAM_PATH,\n",
    "                HMMER_OUTPUT_DIR,\n",
    "                None) for chunk_index,\n",
    "            protein_pair_pid_chunks in enumerate(protein_pair_pid_chunks))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And just like that we ran HMMER against pfam in parallel! The results can be found in the `data` directory. Let's see how the results look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>query_id</th>\n",
       "      <th>accession_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>A0A6L5BYG5</td>\n",
       "      <td>PF03567.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>A0A6L5BTN2</td>\n",
       "      <td>PF00072.27;PF00486.31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A0A6L5BN44</td>\n",
       "      <td>PF00834.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>A0A6L5BWT7</td>\n",
       "      <td>PF02542.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>A0A6L5BWU3</td>\n",
       "      <td>PF04350.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>995</th>\n",
       "      <td>A0A1H3C732</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>A0A1H4LUD9</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>A0A1H3AVF1</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>A0A562J0D1</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>999</th>\n",
       "      <td>A0A1H2ZW83</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       query_id           accession_id\n",
       "0    A0A6L5BYG5             PF03567.17\n",
       "1    A0A6L5BTN2  PF00072.27;PF00486.31\n",
       "2    A0A6L5BN44             PF00834.22\n",
       "3    A0A6L5BWT7             PF02542.19\n",
       "4    A0A6L5BWU3             PF04350.16\n",
       "..          ...                    ...\n",
       "995  A0A1H3C732                    NaN\n",
       "996  A0A1H4LUD9                    NaN\n",
       "997  A0A1H3AVF1                    NaN\n",
       "998  A0A562J0D1                    NaN\n",
       "999  A0A1H2ZW83                    NaN\n",
       "\n",
       "[1000 rows x 2 columns]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.read_csv(f\"{HMMER_OUTPUT_DIR}0_output.csv\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have a function that parses HMMER results! However, we won't show in this notebook, but you can find it in the 'hmmer.py' file in the 'pairpro' package. Please check out the Read the Docs page for more information!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "test_env",
   "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}