analysis/1_fragsys_input_tables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GENERATING MAIN OUTPUT TABLES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains the code used to generate the main output tables from the 1,000s of files that FRAGSYS generates. Due to a size constraint, these files are not in the repository. However, the main tables generated here are, and are the basis of the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## IMPORTING NECESSARY PACKAGES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fragsys_analysis import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## READING INPUT DATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "main_dir = \"/cluster/gjb_lab/2394007/pandda_analysis/phase4/prots\"\n",
    "prot_dirs = sorted(Path(main_dir).iterdir(), key = os.path.getmtime)\n",
    "results_dir = \"./../results/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PROTEIN SUMMARY DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group 1 of P0DTD1 did not present results\n",
      "Group 3 of P0DTD1 did not present results\n"
     ]
    }
   ],
   "source": [
    "results_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\"))) # one for each group\n",
    "    for subdir in subdirs:\n",
    "        res_path = os.path.join(prot_dir, \"{}_{}_results.csv\".format(prot, subdir))\n",
    "        if os.path.isfile(res_path):\n",
    "            prot_res = pd.read_csv(res_path)\n",
    "            results_dfs.append(prot_res)\n",
    "        else:\n",
    "            print(\"Group {} of {} did not present results\".format(subdir, prot))\n",
    "            pass\n",
    "        \n",
    "fragsys_results_df = pd.concat(results_dfs).reset_index(drop = True)\n",
    "fragsys_results_df[\"vars_per_seq\"] = round(fragsys_results_df.n_vars/fragsys_results_df.n_human_seqs, 2)\n",
    "fragsys_results_df[\"vars_per_res\"] = round(fragsys_results_df.n_vars/fragsys_results_df.n_human_res, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(37, 13)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fragsys_results_df.shape # total number of protein segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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>acc</th>\n",
       "      <th>group</th>\n",
       "      <th>n_strucs</th>\n",
       "      <th>n_ligs</th>\n",
       "      <th>n_un_ligs</th>\n",
       "      <th>n_bs</th>\n",
       "      <th>n_seqs</th>\n",
       "      <th>n_human_seqs</th>\n",
       "      <th>n_var_seqs</th>\n",
       "      <th>n_vars</th>\n",
       "      <th>n_human_res</th>\n",
       "      <th>vars_per_seq</th>\n",
       "      <th>vars_per_res</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "      <td>24</td>\n",
       "      <td>25</td>\n",
       "      <td>24</td>\n",
       "      <td>5</td>\n",
       "      <td>108</td>\n",
       "      <td>24</td>\n",
       "      <td>20</td>\n",
       "      <td>1159</td>\n",
       "      <td>2098</td>\n",
       "      <td>48.29</td>\n",
       "      <td>0.55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Q96HY7</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>10</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>185</td>\n",
       "      <td>7</td>\n",
       "      <td>7</td>\n",
       "      <td>2429</td>\n",
       "      <td>4067</td>\n",
       "      <td>347.00</td>\n",
       "      <td>0.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Q9Y2J2</td>\n",
       "      <td>0</td>\n",
       "      <td>51</td>\n",
       "      <td>56</td>\n",
       "      <td>51</td>\n",
       "      <td>10</td>\n",
       "      <td>214</td>\n",
       "      <td>56</td>\n",
       "      <td>54</td>\n",
       "      <td>5860</td>\n",
       "      <td>12756</td>\n",
       "      <td>104.64</td>\n",
       "      <td>0.46</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      acc  group  n_strucs  n_ligs  n_un_ligs  n_bs  n_seqs  n_human_seqs  \\\n",
       "0  P15379      0        24      25         24     5     108            24   \n",
       "1  Q96HY7      0         5      10          5     5     185             7   \n",
       "2  Q9Y2J2      0        51      56         51    10     214            56   \n",
       "\n",
       "   n_var_seqs  n_vars  n_human_res  vars_per_seq  vars_per_res  \n",
       "0          20    1159         2098         48.29          0.55  \n",
       "1           7    2429         4067        347.00          0.60  \n",
       "2          54    5860        12756        104.64          0.46  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fragsys_results_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "fragsys_results_df.to_pickle(os.path.join(results_dir, \"all_prots.pkl\")) # all protein segments in data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BINDING SITES DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group P0DTD1 of 1 did not present binding sites data\n",
      "Group P0DTD1 of 3 did not present binding sites data\n"
     ]
    }
   ],
   "source": [
    "bs_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\")))\n",
    "    for subdir in subdirs:\n",
    "        bs_res_path = os.path.join(prot_dir, \"results\", \"{}_{}_BS_df.csv\".format(prot, subdir))\n",
    "        if os.path.isfile(bs_res_path):\n",
    "            prot_bs_res = pd.read_csv(bs_res_path)\n",
    "            prot_bs_def = pd.read_csv(os.path.join(prot_dir, \"results\", \"{}_{}_BS_def_OC_single_i_rel_0.66.csv\".format(prot, subdir)))\n",
    "            bs_lig_occ_dict = {\"BS\"+str(k):v for k, v in prot_bs_def.binding_site.value_counts().to_dict().items()}\n",
    "            prot_bs_res[\"number_ligs\"] = prot_bs_res.bs_id.map(bs_lig_occ_dict)\n",
    "            prot_bs_res[\"prop_ligs\"] = prot_bs_res.number_ligs/prot_bs_res.number_ligs.sum()\n",
    "            prot_bs_res[\"protein\"] = prot\n",
    "            prot_bs_res[\"group\"] = subdir\n",
    "            bs_dfs.append(prot_bs_res)\n",
    "        else:\n",
    "            print(\"Group {} of {} did not present binding sites data\".format(prot, subdir))\n",
    "            pass\n",
    "all_bss_dfs = round(pd.concat(bs_dfs).reset_index(drop = True), 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(293, 14)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_bss_dfs.shape # total number of residues with structural information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>bs_id</th>\n",
       "      <th>vars</th>\n",
       "      <th>occ</th>\n",
       "      <th>vars_per_occ</th>\n",
       "      <th>MES</th>\n",
       "      <th>p</th>\n",
       "      <th>norm_shenkin_rel</th>\n",
       "      <th>shenkin_ci</th>\n",
       "      <th>MES_ci</th>\n",
       "      <th>number_bs_res</th>\n",
       "      <th>number_ligs</th>\n",
       "      <th>prop_ligs</th>\n",
       "      <th>protein</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BS0</td>\n",
       "      <td>23</td>\n",
       "      <td>37</td>\n",
       "      <td>0.6216</td>\n",
       "      <td>0.1203</td>\n",
       "      <td>0.6838</td>\n",
       "      <td>14.7446</td>\n",
       "      <td>14.3426</td>\n",
       "      <td>0.5255</td>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "      <td>0.08</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BS1</td>\n",
       "      <td>6</td>\n",
       "      <td>14</td>\n",
       "      <td>0.4286</td>\n",
       "      <td>-0.2554</td>\n",
       "      <td>0.8153</td>\n",
       "      <td>7.9722</td>\n",
       "      <td>3.9064</td>\n",
       "      <td>0.9591</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>0.12</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>BS2</td>\n",
       "      <td>108</td>\n",
       "      <td>176</td>\n",
       "      <td>0.6136</td>\n",
       "      <td>0.1153</td>\n",
       "      <td>0.3647</td>\n",
       "      <td>41.1386</td>\n",
       "      <td>14.1057</td>\n",
       "      <td>0.2511</td>\n",
       "      <td>17</td>\n",
       "      <td>18</td>\n",
       "      <td>0.72</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  bs_id  vars  occ  vars_per_occ     MES       p  norm_shenkin_rel  \\\n",
       "0   BS0    23   37        0.6216  0.1203  0.6838           14.7446   \n",
       "1   BS1     6   14        0.4286 -0.2554  0.8153            7.9722   \n",
       "2   BS2   108  176        0.6136  0.1153  0.3647           41.1386   \n",
       "\n",
       "   shenkin_ci  MES_ci  number_bs_res  number_ligs  prop_ligs protein group  \n",
       "0     14.3426  0.5255              8            2       0.08  P15379     0  \n",
       "1      3.9064  0.9591              5            3       0.12  P15379     0  \n",
       "2     14.1057  0.2511             17           18       0.72  P15379     0  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_bss_dfs.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "all_bss_dfs.to_pickle(os.path.join(results_dir, \"all_bss.pkl\")) # all binding sites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BINDING RESIDUES DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group P0DTD1 of 1 did not present results\n",
      "Group P0DTD1 of 3 did not present results\n"
     ]
    }
   ],
   "source": [
    "fragsys_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\")))\n",
    "    for subdir in subdirs:\n",
    "        res_path = os.path.join(prot_dir, \"results\", \"{}_{}_fragsys_df.csv\".format(prot, subdir))\n",
    "        if os.path.isfile(res_path):\n",
    "            prot_res = pd.read_csv(res_path)\n",
    "            prot_res[\"protein\"] = prot\n",
    "            prot_res[\"group\"] = subdir\n",
    "            fragsys_dfs.append(prot_res)\n",
    "        else:\n",
    "            print(\"Group {} of {} did not present results\".format(prot, subdir))\n",
    "            pass\n",
    "        \n",
    "fragsyss_df = pd.concat(fragsys_dfs).reset_index(drop = True)\n",
    "fragsyss_df.SS = fragsyss_df.SS.fillna(\"C\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14172, 62)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fragsyss_df.shape # total number of ligand binding residues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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>index</th>\n",
       "      <th>UniProt_ResNum</th>\n",
       "      <th>Pfam_column</th>\n",
       "      <th>alignment_column</th>\n",
       "      <th>shenkin</th>\n",
       "      <th>occ</th>\n",
       "      <th>gaps</th>\n",
       "      <th>occ_pct</th>\n",
       "      <th>gaps_pct</th>\n",
       "      <th>rel_norm_shenkin</th>\n",
       "      <th>...</th>\n",
       "      <th>BS14</th>\n",
       "      <th>BS15</th>\n",
       "      <th>BS16</th>\n",
       "      <th>BS17</th>\n",
       "      <th>BS18</th>\n",
       "      <th>BS19</th>\n",
       "      <th>BS20</th>\n",
       "      <th>BS21</th>\n",
       "      <th>BS22</th>\n",
       "      <th>BS23</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8</td>\n",
       "      <td>31</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>51.846631</td>\n",
       "      <td>57</td>\n",
       "      <td>51</td>\n",
       "      <td>0.527778</td>\n",
       "      <td>0.472222</td>\n",
       "      <td>82.865612</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13</td>\n",
       "      <td>36</td>\n",
       "      <td>18</td>\n",
       "      <td>18</td>\n",
       "      <td>14.127808</td>\n",
       "      <td>97</td>\n",
       "      <td>11</td>\n",
       "      <td>0.898148</td>\n",
       "      <td>0.101852</td>\n",
       "      <td>14.690627</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>35</td>\n",
       "      <td>58</td>\n",
       "      <td>44</td>\n",
       "      <td>44</td>\n",
       "      <td>58.025383</td>\n",
       "      <td>107</td>\n",
       "      <td>1</td>\n",
       "      <td>0.990741</td>\n",
       "      <td>0.009259</td>\n",
       "      <td>94.033413</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 62 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   index  UniProt_ResNum  Pfam_column  alignment_column    shenkin  occ  gaps  \\\n",
       "0      8              31           11                11  51.846631   57    51   \n",
       "1     13              36           18                18  14.127808   97    11   \n",
       "2     35              58           44                44  58.025383  107     1   \n",
       "\n",
       "    occ_pct  gaps_pct  rel_norm_shenkin  ...  BS14  BS15  BS16  BS17  BS18  \\\n",
       "0  0.527778  0.472222         82.865612  ...   NaN   NaN   NaN   NaN   NaN   \n",
       "1  0.898148  0.101852         14.690627  ...   NaN   NaN   NaN   NaN   NaN   \n",
       "2  0.990741  0.009259         94.033413  ...   NaN   NaN   NaN   NaN   NaN   \n",
       "\n",
       "   BS19  BS20  BS21  BS22  BS23  \n",
       "0   NaN   NaN   NaN   NaN   NaN  \n",
       "1   NaN   NaN   NaN   NaN   NaN  \n",
       "2   NaN   NaN   NaN   NaN   NaN  \n",
       "\n",
       "[3 rows x 62 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fragsyss_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fragsyss_df.to_pickle(os.path.join(results_dir, \"all_bs_ress.pkl\")) # all binding site residues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DSSP DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "dssp_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\")))\n",
    "    for subdir in subdirs:\n",
    "        dssp_subdir = os.path.join(prot_dir, \"results\", \"dssp\", subdir)\n",
    "        sifts_subdir = os.path.join(prot_dir, \"results\", \"sifts\", subdir)\n",
    "        dssp_files = [file for file in os.listdir(dssp_subdir) if file.startswith(\"dssp\")]\n",
    "        pdb_id = dssp_files[0].split(\"_\")[1]\n",
    "        example_dssp = os.path.join(dssp_subdir, dssp_files[0]) #only grabbing one example dssp file\n",
    "        example_sifts = os.path.join(sifts_subdir, \"sifts_mapping_{}.csv\".format(pdb_id))\n",
    "        if os.path.isfile(example_dssp):\n",
    "            dssp_df = pd.read_csv(example_dssp)\n",
    "            sifts_df = pd.read_csv(example_sifts)\n",
    "            dssp_df[\"protein\"] = prot\n",
    "            dssp_df[\"group\"] = subdir\n",
    "            dssp_df[\"structure\"] = pdb_id\n",
    "            merged_df = dssp_df.merge(sifts_df[[\"PDB_ResNum\", \"PDB_ChainID\", \"UniProt_ResNum\"]], how = \"left\", left_on = [\"PDB_ResNum\", \"CHAIN\"], right_on = [\"PDB_ResNum\", \"PDB_ChainID\"])\n",
    "            dssp_dfs.append(merged_df)\n",
    "        else:\n",
    "            print(\"Group {} of {} did not present DSSP data\".format(prot, subdir))\n",
    "            pass\n",
    "        \n",
    "all_dssp_dfs = pd.concat(dssp_dfs)\n",
    "all_dssp_dfs_filt = all_dssp_dfs.copy().query('UniProt_ResNum == UniProt_ResNum')\n",
    "all_dssp_dfs_filt.UniProt_ResNum = all_dssp_dfs_filt.UniProt_ResNum.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(13016, 17)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_dssp_dfs_filt.shape # total number of residues with structural information"
   ]
  },
  {
   "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>PDB_ResNum</th>\n",
       "      <th>CHAIN</th>\n",
       "      <th>AA</th>\n",
       "      <th>SS</th>\n",
       "      <th>ACC</th>\n",
       "      <th>TCO</th>\n",
       "      <th>KAPPA</th>\n",
       "      <th>ALPHA</th>\n",
       "      <th>PHI</th>\n",
       "      <th>PSI</th>\n",
       "      <th>CHAIN_FULL</th>\n",
       "      <th>RSA</th>\n",
       "      <th>protein</th>\n",
       "      <th>group</th>\n",
       "      <th>structure</th>\n",
       "      <th>PDB_ChainID</th>\n",
       "      <th>UniProt_ResNum</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>25</td>\n",
       "      <td>A</td>\n",
       "      <td>Q</td>\n",
       "      <td>E</td>\n",
       "      <td>82</td>\n",
       "      <td>-0.950</td>\n",
       "      <td>360.0</td>\n",
       "      <td>-166.5</td>\n",
       "      <td>-135.2</td>\n",
       "      <td>146.3</td>\n",
       "      <td>A</td>\n",
       "      <td>41.414</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "      <td>5sc6</td>\n",
       "      <td>A</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>26</td>\n",
       "      <td>A</td>\n",
       "      <td>I</td>\n",
       "      <td>E</td>\n",
       "      <td>4</td>\n",
       "      <td>-0.999</td>\n",
       "      <td>10.9</td>\n",
       "      <td>-155.9</td>\n",
       "      <td>-129.4</td>\n",
       "      <td>128.5</td>\n",
       "      <td>A</td>\n",
       "      <td>2.367</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "      <td>5sc6</td>\n",
       "      <td>A</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>27</td>\n",
       "      <td>A</td>\n",
       "      <td>D</td>\n",
       "      <td>E</td>\n",
       "      <td>37</td>\n",
       "      <td>-0.900</td>\n",
       "      <td>6.8</td>\n",
       "      <td>-167.5</td>\n",
       "      <td>-101.2</td>\n",
       "      <td>132.1</td>\n",
       "      <td>A</td>\n",
       "      <td>22.699</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "      <td>5sc6</td>\n",
       "      <td>A</td>\n",
       "      <td>25</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   PDB_ResNum CHAIN AA SS  ACC    TCO  KAPPA  ALPHA    PHI    PSI CHAIN_FULL  \\\n",
       "1          25     A  Q  E   82 -0.950  360.0 -166.5 -135.2  146.3          A   \n",
       "2          26     A  I  E    4 -0.999   10.9 -155.9 -129.4  128.5          A   \n",
       "3          27     A  D  E   37 -0.900    6.8 -167.5 -101.2  132.1          A   \n",
       "\n",
       "      RSA protein group structure PDB_ChainID  UniProt_ResNum  \n",
       "1  41.414  P15379     0      5sc6           A              23  \n",
       "2   2.367  P15379     0      5sc6           A              24  \n",
       "3  22.699  P15379     0      5sc6           A              25  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_dssp_dfs_filt.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_dssp_dfs_filt.to_pickle(os.path.join(results_dir, \"all_dssp_dfs.pkl\")) # all dssp data frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## VARIANTS DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group P0DTD1 of 1 did not present variation data\n",
      "Group P0DTD1 of 3 did not present variation data\n"
     ]
    }
   ],
   "source": [
    "missense_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\")))\n",
    "    for subdir in subdirs:\n",
    "        missense_df_path = os.path.join(prot_dir,\"results/varalign/{}/{}_{}_missense_df.csv\".format(subdir, prot, subdir))\n",
    "        if os.path.isfile(missense_df_path):\n",
    "            miss_df = pd.read_csv(missense_df_path)\n",
    "            miss_df[\"protein\"] = prot\n",
    "            miss_df[\"group\"] = subdir\n",
    "            missense_dfs.append(miss_df)\n",
    "        else:\n",
    "            print(\"Group {} of {} did not present variation data\".format(prot, subdir))\n",
    "            pass\n",
    "        \n",
    "all_missense_df = pd.concat(missense_dfs).reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10408, 17)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_missense_df.shape # total number of alignment columns with variation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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>col</th>\n",
       "      <th>shenkin</th>\n",
       "      <th>occ</th>\n",
       "      <th>gaps</th>\n",
       "      <th>occ_pct</th>\n",
       "      <th>gaps_pct</th>\n",
       "      <th>variants</th>\n",
       "      <th>rel_norm_shenkin</th>\n",
       "      <th>abs_norm_shenkin</th>\n",
       "      <th>oddsratio</th>\n",
       "      <th>log_oddsratio</th>\n",
       "      <th>pvalue</th>\n",
       "      <th>ci_dist</th>\n",
       "      <th>miss_class</th>\n",
       "      <th>miss_color</th>\n",
       "      <th>protein</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>1</td>\n",
       "      <td>20</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.810181</td>\n",
       "      <td>0.593427</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2.772787</td>\n",
       "      <td>CME</td>\n",
       "      <td>green</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>9.639785</td>\n",
       "      <td>3</td>\n",
       "      <td>18</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.90</td>\n",
       "      <td>2</td>\n",
       "      <td>6.578738</td>\n",
       "      <td>3.192794</td>\n",
       "      <td>1.206678</td>\n",
       "      <td>0.187871</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.790666</td>\n",
       "      <td>CME</td>\n",
       "      <td>green</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>11.927558</td>\n",
       "      <td>4</td>\n",
       "      <td>17</td>\n",
       "      <td>0.15</td>\n",
       "      <td>0.85</td>\n",
       "      <td>3</td>\n",
       "      <td>10.713779</td>\n",
       "      <td>5.199612</td>\n",
       "      <td>1.358038</td>\n",
       "      <td>0.306041</td>\n",
       "      <td>0.704905</td>\n",
       "      <td>1.498695</td>\n",
       "      <td>CME</td>\n",
       "      <td>green</td>\n",
       "      <td>P15379</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   col    shenkin  occ  gaps  occ_pct  gaps_pct  variants  rel_norm_shenkin  \\\n",
       "0    1   6.000000    1    20     0.00      1.00         1          0.000000   \n",
       "1    2   9.639785    3    18     0.10      0.90         2          6.578738   \n",
       "2    3  11.927558    4    17     0.15      0.85         3         10.713779   \n",
       "\n",
       "   abs_norm_shenkin  oddsratio  log_oddsratio    pvalue   ci_dist miss_class  \\\n",
       "0          0.000000   1.810181       0.593427  1.000000  2.772787        CME   \n",
       "1          3.192794   1.206678       0.187871  1.000000  1.790666        CME   \n",
       "2          5.199612   1.358038       0.306041  0.704905  1.498695        CME   \n",
       "\n",
       "  miss_color protein group  \n",
       "0      green  P15379     0  \n",
       "1      green  P15379     0  \n",
       "2      green  P15379     0  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_missense_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_missense_df.to_pickle(os.path.join(results_dir, \"all_miss_dfs.pkl\")) # all variant-containing columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## INTERACTIONS DATAFRAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cons_dfs = []\n",
    "for prot_dir in prot_dirs:\n",
    "    prot = str(prot_dir).split(\"/\")[-1]\n",
    "    subdirs = sorted(os.listdir(os.path.join(prot_dir, \"unsupp_cifs\")))\n",
    "    for subdir in subdirs:\n",
    "        arpeggio_subdir = os.path.join(prot_dir, \"results\", \"arpeggio\", subdir)\n",
    "        arpeggio_files = [file for file in os.listdir(arpeggio_subdir) if file.startswith(\"arpeggio_all_cons_split\")]\n",
    "        for file in arpeggio_files:\n",
    "            pdb_id = file[-8:-4]\n",
    "            file_df = pd.read_csv(os.path.join(arpeggio_subdir, file))\n",
    "            file_df[\"struc\"] = pdb_id\n",
    "            file_df[\"protein\"] = prot\n",
    "            all_cons_dfs.append(file_df)\n",
    "all_fragsys_cons = pd.concat(all_cons_dfs).reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(269013, 30)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_fragsys_cons.shape # total number of interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>Chain (Atom1)</th>\n",
       "      <th>ResNum (Atom1)</th>\n",
       "      <th>ResName (Atom1)</th>\n",
       "      <th>Atom (Atom1)</th>\n",
       "      <th>Chain (Atom2)</th>\n",
       "      <th>ResNum (Atom2)</th>\n",
       "      <th>ResName (Atom2)</th>\n",
       "      <th>Atom (Atom2)</th>\n",
       "      <th>Clash</th>\n",
       "      <th>Covalent</th>\n",
       "      <th>...</th>\n",
       "      <th>Carbonyl</th>\n",
       "      <th>Polar</th>\n",
       "      <th>Weak Polar</th>\n",
       "      <th>Atom proximity</th>\n",
       "      <th>Vdw proximity</th>\n",
       "      <th>Interacting entities</th>\n",
       "      <th>contact_type</th>\n",
       "      <th>UniProt_Resnum</th>\n",
       "      <th>struc</th>\n",
       "      <th>protein</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>A</td>\n",
       "      <td>33</td>\n",
       "      <td>ARG</td>\n",
       "      <td>CD</td>\n",
       "      <td>B</td>\n",
       "      <td>201</td>\n",
       "      <td>NW4</td>\n",
       "      <td>O2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3.771</td>\n",
       "      <td>0.551</td>\n",
       "      <td>INTER</td>\n",
       "      <td>sidechain</td>\n",
       "      <td>31.0</td>\n",
       "      <td>5sbv</td>\n",
       "      <td>P15379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>A</td>\n",
       "      <td>33</td>\n",
       "      <td>ARG</td>\n",
       "      <td>CZ</td>\n",
       "      <td>B</td>\n",
       "      <td>201</td>\n",
       "      <td>NW4</td>\n",
       "      <td>O2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4.749</td>\n",
       "      <td>1.529</td>\n",
       "      <td>INTER</td>\n",
       "      <td>sidechain</td>\n",
       "      <td>31.0</td>\n",
       "      <td>5sbv</td>\n",
       "      <td>P15379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A</td>\n",
       "      <td>33</td>\n",
       "      <td>ARG</td>\n",
       "      <td>NE</td>\n",
       "      <td>B</td>\n",
       "      <td>201</td>\n",
       "      <td>NW4</td>\n",
       "      <td>O2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3.821</td>\n",
       "      <td>0.751</td>\n",
       "      <td>INTER</td>\n",
       "      <td>sidechain</td>\n",
       "      <td>31.0</td>\n",
       "      <td>5sbv</td>\n",
       "      <td>P15379</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 30 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "  Chain (Atom1)  ResNum (Atom1) ResName (Atom1) Atom (Atom1) Chain (Atom2)  \\\n",
       "0             A              33             ARG           CD             B   \n",
       "1             A              33             ARG           CZ             B   \n",
       "2             A              33             ARG           NE             B   \n",
       "\n",
       "   ResNum (Atom2) ResName (Atom2) Atom (Atom2)  Clash  Covalent  ...  \\\n",
       "0             201             NW4           O2      0         0  ...   \n",
       "1             201             NW4           O2      0         0  ...   \n",
       "2             201             NW4           O2      0         0  ...   \n",
       "\n",
       "   Carbonyl  Polar  Weak Polar  Atom proximity  Vdw proximity  \\\n",
       "0         0      0           0           3.771          0.551   \n",
       "1         0      0           0           4.749          1.529   \n",
       "2         0      0           0           3.821          0.751   \n",
       "\n",
       "   Interacting entities  contact_type  UniProt_Resnum  struc  protein  \n",
       "0                 INTER     sidechain            31.0   5sbv   P15379  \n",
       "1                 INTER     sidechain            31.0   5sbv   P15379  \n",
       "2                 INTER     sidechain            31.0   5sbv   P15379  \n",
       "\n",
       "[3 rows x 30 columns]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_fragsys_cons.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_fragsys_cons.to_pickle(os.path.join(results_dir, \"all_bss_cons.pkl\")) # all binding site contacts"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:myenv2]",
   "language": "python",
   "name": "conda-env-myenv2-py"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}