notebooks/1-CreateDataset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Dataset\n",
    "This notebook extracts from the Protein Data Bank information about the secondary structure of proteins. The ultimate goal is to assign a fold classification from a protein sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**Rule 2: Document the Process, Not Just the Results.** Here we describe the steps how to produce the dataset.\n",
    "\n",
    "\n",
    "**Rule 7: Build a Pipeline.** Besides documenting all steps, the entire process of dataset creation from the original data files in the /data directory is automated. There are no manual steps.\n",
    "\n",
    "\n",
    "**Rule 8: Share and Explain Your Data.** To enable reproducibility we provide a /data directory with data files and a file that describes the datasets with download locations and dates.\n",
    "\n",
    "\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# column names\n",
    "value_col = \"foldClass\" # fold class to be predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np     \n",
    "import pdbutils "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read a Representative Set of Protein Chains\n",
    "Protein sequences in the Protein Data Bank (PDB) are redundant. For example, there are more than 300 structures of hemoglobin in the PDB. To reduce biases in our analysis, we use a representative set of protein chains, i.e., a set of protein chains with minimal sequence identity among its members. For this analysis we downloaded a non-redundant set from the [PISCES website](http://dunbrack.fccc.edu/PISCES.php) with a maximum of 25% sequence identity and an X-ray resolution of <= 3.0 &#197;. \n",
    "\n",
    "Wang G, Dunbrack RL Jr (2005) PISCES: recent improvements to a PDB sequence culling server, Nucleic Acids Res. 33, W94-8. [doi: 10.1093/nar/gki402](https://doi.org/10.1093/nar/gki402)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of representative PDB chains: 14051 \n",
      "\n"
     ]
    },
    {
     "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>length</th>\n",
       "      <th>Exptl.</th>\n",
       "      <th>resolution</th>\n",
       "      <th>R-factor</th>\n",
       "      <th>FreeRvalue</th>\n",
       "      <th>pdbChainId</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>330</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.29</td>\n",
       "      <td>12AS.A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>366</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.1</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.26</td>\n",
       "      <td>16VP.A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>348</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.6</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.34</td>\n",
       "      <td>1A0I.A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>413</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>1.7</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.22</td>\n",
       "      <td>1A12.C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>108</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.25</td>\n",
       "      <td>1A1X.A</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   length Exptl.  resolution  R-factor  FreeRvalue pdbChainId\n",
       "0     330   XRAY         2.2      0.16        0.29     12AS.A\n",
       "1     366   XRAY         2.1      0.19        0.26     16VP.A\n",
       "2     348   XRAY         2.6      0.22        0.34     1A0I.A\n",
       "3     413   XRAY         1.7      0.19        0.22     1A12.C\n",
       "4     108   XRAY         2.0      0.21        0.25     1A1X.A"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "representatives = pdbutils.read_pisces_representatives('../data/cullpdb_pc25_res3.0_R1.0_d180920_chains14051.gz')\n",
    "print(\"Number of representative PDB chains:\", representatives.shape[0], \"\\n\")\n",
    "representatives.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read Protein Sequence and Secondary Structure Data\n",
    "Secondary structure of proteins is most commonly assigned with the DSSP method.\n",
    "\n",
    "Kabsch W, Sander C (1983) Dictionary of protein secondary structure: Pattern recognition of hydrogen‐bonded and geometrical features. Biopolymers 22, 2577–637. [doi 0.1002/bip.360221211](https://doi.org/10.1002/bip.360221211)\n",
    "\n",
    "DSSP defines [8 classes of secondary structure](https://en.wikipedia.org/wiki/Protein_secondary_structure):\n",
    "\n",
    "| DSSP Class  | Code |\n",
    "| :---        | ---: |\n",
    "| 4-turn helix (&#945; helix), min. length 4 residues | H |\n",
    "| Isolated &#946;-bridge (single pair &#946;-sheet hydrogen bond formation) | B |\n",
    "| Extended strand in parallel or anti-parallel &#946;-sheet conformation, min length 2 residues | E |\n",
    "| 3-turn helix (3<sub>10</sub> helix), min. length 3 residues | G |\n",
    "| 5-turn helix (&#960; helix), min. length 5 residues | I |\n",
    "| Hydrogen bonded turn (3, 4 or 5 turn) | T |\n",
    "| Bend (the only non-hydrogen-bond based assignment) | S |\n",
    "| Coil (residues which are not in any of the above conformations) | C |\n",
    "\n",
    "The [RCSB Protein Data Bank](https://www.rcsb.org) provides protein sequences and DSSP secondary structure assignments. The method below reads a local copy of this file and returns the data as a Pandas dataframe. \n",
    "\n",
    "The **secondary_structure** string below is the secondary structure assignment for each amino acid residue in the **sequence**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of protein chains in PDB: 402007 \n",
      "\n"
     ]
    },
    {
     "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>pdbChainId</th>\n",
       "      <th>sequence</th>\n",
       "      <th>secondary_structure</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>101M.A</td>\n",
       "      <td>MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>102L.A</td>\n",
       "      <td>MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE...</td>\n",
       "      <td>CCHHHHHHHHHCCEEEEEECTTSCEEEETTEEEESSSCTTTHHHHH...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>102M.A</td>\n",
       "      <td>MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>103L.A</td>\n",
       "      <td>MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNSLDAAK...</td>\n",
       "      <td>CCHHHHHHHHHCCEEEEEECTTSCEEEETTEECCCCCCCCCHHHHH...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>103M.A</td>\n",
       "      <td>MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pdbChainId                                           sequence  \\\n",
       "0     101M.A  MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...   \n",
       "1     102L.A  MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE...   \n",
       "2     102M.A  MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...   \n",
       "3     103L.A  MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNSLDAAK...   \n",
       "4     103M.A  MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...   \n",
       "\n",
       "                                 secondary_structure  \n",
       "0  CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...  \n",
       "1  CCHHHHHHHHHCCEEEEEECTTSCEEEETTEEEESSSCTTTHHHHH...  \n",
       "2  CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...  \n",
       "3  CCHHHHHHHHHCCEEEEEECTTSCEEEETTEECCCCCCCCCHHHHH...  \n",
       "4  CCCCHHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHCGGGGGGCTT...  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sec_struct = pdbutils.read_secondary_structure('../data/ss_dis.txt.gz')\n",
    "print(\"Number of protein chains in PDB:\", sec_struct.shape[0], \"\\n\")\n",
    "sec_struct.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find the Intersection between the two Datasets\n",
    "Note the high level of redundancy in the PDB: the representative set contains about 14,000 protein chains, whereas the entire PDB contains more than 400,000 protein chains.\n",
    "\n",
    "By merging the representative set of protein chains with the secondary structure dataset we obtain the intersection between the two sets. Both datasets contain a unique identifier column **pdbChainId** that is used to join the datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of representative chains with secondary structure data: 13791\n"
     ]
    },
    {
     "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>pdbChainId</th>\n",
       "      <th>sequence</th>\n",
       "      <th>secondary_structure</th>\n",
       "      <th>length</th>\n",
       "      <th>Exptl.</th>\n",
       "      <th>resolution</th>\n",
       "      <th>R-factor</th>\n",
       "      <th>FreeRvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>12AS.A</td>\n",
       "      <td>MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQD...</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHHHHHHHHHHCEEECCCCSEEETTSSCSC...</td>\n",
       "      <td>330</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>16VP.A</td>\n",
       "      <td>SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...</td>\n",
       "      <td>CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...</td>\n",
       "      <td>366</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.1</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1A0I.A</td>\n",
       "      <td>VNIKTNPFKAVSFVESAIKKALDNAGYLIAEIKYDGVRGNICVDNT...</td>\n",
       "      <td>CTTCCCCEEEEECCHHHHHHHHHHHSSEEEEECCCSEEEEEEEETT...</td>\n",
       "      <td>348</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.6</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1A12.C</td>\n",
       "      <td>RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...</td>\n",
       "      <td>CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...</td>\n",
       "      <td>413</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>1.7</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1A1X.A</td>\n",
       "      <td>GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...</td>\n",
       "      <td>CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...</td>\n",
       "      <td>108</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pdbChainId                                           sequence  \\\n",
       "0     12AS.A  MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQD...   \n",
       "1     16VP.A  SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...   \n",
       "2     1A0I.A  VNIKTNPFKAVSFVESAIKKALDNAGYLIAEIKYDGVRGNICVDNT...   \n",
       "3     1A12.C  RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...   \n",
       "4     1A1X.A  GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...   \n",
       "\n",
       "                                 secondary_structure  length Exptl.  \\\n",
       "0  CCCCHHHHHHHHHHHHHHHHHHHHHHHCEEECCCCSEEETTSSCSC...     330   XRAY   \n",
       "1  CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...     366   XRAY   \n",
       "2  CTTCCCCEEEEECCHHHHHHHHHHHSSEEEEECCCSEEEEEEEETT...     348   XRAY   \n",
       "3  CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...     413   XRAY   \n",
       "4  CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...     108   XRAY   \n",
       "\n",
       "   resolution  R-factor  FreeRvalue  \n",
       "0         2.2      0.16        0.29  \n",
       "1         2.1      0.19        0.26  \n",
       "2         2.6      0.22        0.34  \n",
       "3         1.7      0.19        0.22  \n",
       "4         2.0      0.21        0.25  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = sec_struct.merge(representatives, left_on='pdbChainId', right_on='pdbChainId', how='inner')\n",
    "print(\"Number of representative chains with secondary structure data:\", df.shape[0])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Secondary Structure Content\n",
    "To reduce the DSSP 8-state classification to a 3-state classification we group related secondary structure elements: alpha (I, H, G), beta (E, B), and coil (S, T, C). Then we calculate the fraction of the amino acid residues in each of the three classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>pdbChainId</th>\n",
       "      <th>sequence</th>\n",
       "      <th>secondary_structure</th>\n",
       "      <th>length</th>\n",
       "      <th>Exptl.</th>\n",
       "      <th>resolution</th>\n",
       "      <th>R-factor</th>\n",
       "      <th>FreeRvalue</th>\n",
       "      <th>alpha</th>\n",
       "      <th>beta</th>\n",
       "      <th>coil</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>12AS.A</td>\n",
       "      <td>MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQD...</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHHHHHHHHHHCEEECCCCSEEETTSSCSC...</td>\n",
       "      <td>330</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.29</td>\n",
       "      <td>0.345455</td>\n",
       "      <td>0.206061</td>\n",
       "      <td>0.448485</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>16VP.A</td>\n",
       "      <td>SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...</td>\n",
       "      <td>CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...</td>\n",
       "      <td>366</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.1</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.469945</td>\n",
       "      <td>0.046448</td>\n",
       "      <td>0.483607</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1A0I.A</td>\n",
       "      <td>VNIKTNPFKAVSFVESAIKKALDNAGYLIAEIKYDGVRGNICVDNT...</td>\n",
       "      <td>CTTCCCCEEEEECCHHHHHHHHHHHSSEEEEECCCSEEEEEEEETT...</td>\n",
       "      <td>348</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.6</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.34</td>\n",
       "      <td>0.232759</td>\n",
       "      <td>0.318966</td>\n",
       "      <td>0.448276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1A12.C</td>\n",
       "      <td>RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...</td>\n",
       "      <td>CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...</td>\n",
       "      <td>413</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>1.7</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.038741</td>\n",
       "      <td>0.418886</td>\n",
       "      <td>0.542373</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1A1X.A</td>\n",
       "      <td>GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...</td>\n",
       "      <td>CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...</td>\n",
       "      <td>108</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.037037</td>\n",
       "      <td>0.472222</td>\n",
       "      <td>0.490741</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pdbChainId                                           sequence  \\\n",
       "0     12AS.A  MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQD...   \n",
       "1     16VP.A  SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...   \n",
       "2     1A0I.A  VNIKTNPFKAVSFVESAIKKALDNAGYLIAEIKYDGVRGNICVDNT...   \n",
       "3     1A12.C  RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...   \n",
       "4     1A1X.A  GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...   \n",
       "\n",
       "                                 secondary_structure  length Exptl.  \\\n",
       "0  CCCCHHHHHHHHHHHHHHHHHHHHHHHCEEECCCCSEEETTSSCSC...     330   XRAY   \n",
       "1  CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...     366   XRAY   \n",
       "2  CTTCCCCEEEEECCHHHHHHHHHHHSSEEEEECCCSEEEEEEEETT...     348   XRAY   \n",
       "3  CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...     413   XRAY   \n",
       "4  CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...     108   XRAY   \n",
       "\n",
       "   resolution  R-factor  FreeRvalue     alpha      beta      coil  \n",
       "0         2.2      0.16        0.29  0.345455  0.206061  0.448485  \n",
       "1         2.1      0.19        0.26  0.469945  0.046448  0.483607  \n",
       "2         2.6      0.22        0.34  0.232759  0.318966  0.448276  \n",
       "3         1.7      0.19        0.22  0.038741  0.418886  0.542373  \n",
       "4         2.0      0.21        0.25  0.037037  0.472222  0.490741  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def alpha_fraction(s):\n",
    "    return (s.count('I') + s.count('H') + s.count('G')) / len(s)\n",
    "\n",
    "\n",
    "def beta_fraction(s):\n",
    "    return (s.count('E') + s.count('B')) / len(s)\n",
    "\n",
    "\n",
    "def coil_fraction(s):\n",
    "    return (s.count('S') + s.count('T') + s.count('C')) / len(s)\n",
    "\n",
    "\n",
    "df['alpha'] = df.secondary_structure.apply(alpha_fraction)                                                                               \n",
    "df['beta'] = df.secondary_structure.apply(beta_fraction)                                                                           \n",
    "df['coil'] = df.secondary_structure.apply(coil_fraction)                                                                 \n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classify Sequences by Secondary Structure Content\n",
    "Next we classify each protein chain into one of four classes. We use a threshold of 25% to define a predominant class.\n",
    "\n",
    "* alpha: predominantly alpha (> 25%)\n",
    "* beta: predominantly beta (> 25%)\n",
    "* alpha+beta: significant alpha (> 25%) and beta (> 25%)\n",
    "* other: cases that do not fit into the 3 classes above \n",
    "\n",
    "Protein chains in the **other** class will be ignored in the subsequent analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def protein_fold_class(data, min_threshold, max_threshold):\n",
    "    \"\"\"\n",
    "    Returns fold classification:\n",
    "    \"alpha\", \"beta\", \"alpha+beta\", and \"other\" based upon the \n",
    "    fraction of alpha and beta.\n",
    "    \"\"\"\n",
    "    if data.alpha > max_threshold and data.beta < min_threshold:                                \n",
    "        return \"alpha\"                                                                        \n",
    "    elif data.beta > max_threshold and data.alpha < min_threshold:                              \n",
    "        return \"beta\"                                                                         \n",
    "    elif data.alpha > max_threshold and data.beta > max_threshold:                              \n",
    "        return \"alpha+beta\"                                                                   \n",
    "    else:                                                                                     \n",
    "        return \"other\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset size 5370\n"
     ]
    },
    {
     "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>pdbChainId</th>\n",
       "      <th>sequence</th>\n",
       "      <th>secondary_structure</th>\n",
       "      <th>length</th>\n",
       "      <th>Exptl.</th>\n",
       "      <th>resolution</th>\n",
       "      <th>R-factor</th>\n",
       "      <th>FreeRvalue</th>\n",
       "      <th>alpha</th>\n",
       "      <th>beta</th>\n",
       "      <th>coil</th>\n",
       "      <th>foldClass</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>16VP.A</td>\n",
       "      <td>SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...</td>\n",
       "      <td>CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...</td>\n",
       "      <td>366</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.10</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.469945</td>\n",
       "      <td>0.046448</td>\n",
       "      <td>0.483607</td>\n",
       "      <td>alpha</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1A12.C</td>\n",
       "      <td>RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...</td>\n",
       "      <td>CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...</td>\n",
       "      <td>413</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>1.70</td>\n",
       "      <td>0.19</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.038741</td>\n",
       "      <td>0.418886</td>\n",
       "      <td>0.542373</td>\n",
       "      <td>beta</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1A1X.A</td>\n",
       "      <td>GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...</td>\n",
       "      <td>CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...</td>\n",
       "      <td>108</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.00</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.037037</td>\n",
       "      <td>0.472222</td>\n",
       "      <td>0.490741</td>\n",
       "      <td>beta</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1A2X.B</td>\n",
       "      <td>GDEEKRNRAITARRQHLKSVMLQIAATELEKEEGRREAEKQNYLAEH</td>\n",
       "      <td>CCCCHHHHHHHHHHHHHHHHHHHHHHHHHHHTCCCCCCCCCCCCCCC</td>\n",
       "      <td>47</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>2.30</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.33</td>\n",
       "      <td>0.574468</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.425532</td>\n",
       "      <td>alpha</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1A62.A</td>\n",
       "      <td>MNLTELKNTPVSELITLGENMGLENLARMRKQDIIFAILKQHAKSG...</td>\n",
       "      <td>CBHHHHHTSCHHHHHHHHHTTTCCCCTTSCHHHHHHHHHHHHHHTT...</td>\n",
       "      <td>130</td>\n",
       "      <td>XRAY</td>\n",
       "      <td>1.55</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.284615</td>\n",
       "      <td>0.261538</td>\n",
       "      <td>0.453846</td>\n",
       "      <td>alpha+beta</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  pdbChainId                                           sequence  \\\n",
       "1     16VP.A  SRMPSPPMPVPPAALFNRLLDDLGFSAGPALCTMLDTWNEDLFSAL...   \n",
       "3     1A12.C  RRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...   \n",
       "4     1A1X.A  GSAGEDVGAPPDHLWVHQEGIYRDEYQRTWVAVVEEETSFLRARVQ...   \n",
       "5     1A2X.B    GDEEKRNRAITARRQHLKSVMLQIAATELEKEEGRREAEKQNYLAEH   \n",
       "8     1A62.A  MNLTELKNTPVSELITLGENMGLENLARMRKQDIIFAILKQHAKSG...   \n",
       "\n",
       "                                 secondary_structure  length Exptl.  \\\n",
       "1  CCSCCCCCCCCHHHHHHHHHHHHTCTTHHHHHHHHHHCCCCCSTTS...     366   XRAY   \n",
       "3  CCCCCCCCCCCCCCCCCCCTTCCCCCBEEEEEEECTTSTTCSCTTC...     413   XRAY   \n",
       "4  CCCCCCCCCCCSEEEEEETTEEEETTSCEEEEEEEECSSCEEEEEE...     108   XRAY   \n",
       "5    CCCCHHHHHHHHHHHHHHHHHHHHHHHHHHHTCCCCCCCCCCCCCCC      47   XRAY   \n",
       "8  CBHHHHHTSCHHHHHHHHHTTTCCCCTTSCHHHHHHHHHHHHHHTT...     130   XRAY   \n",
       "\n",
       "   resolution  R-factor  FreeRvalue     alpha      beta      coil   foldClass  \n",
       "1        2.10      0.19        0.26  0.469945  0.046448  0.483607       alpha  \n",
       "3        1.70      0.19        0.22  0.038741  0.418886  0.542373        beta  \n",
       "4        2.00      0.21        0.25  0.037037  0.472222  0.490741        beta  \n",
       "5        2.30      0.22        0.33  0.574468  0.000000  0.425532       alpha  \n",
       "8        1.55      0.22        0.25  0.284615  0.261538  0.453846  alpha+beta  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# assign protein fold class\n",
    "df[value_col] = df.apply(protein_fold_class, min_threshold=0.05, max_threshold=0.25, axis=1)\n",
    "\n",
    "# exclude protein chains without a dominant classification from further analysis.\n",
    "df = df[df[value_col] != 'other']\n",
    "\n",
    "print(\"Dataset size\", df.shape[0])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save Dataset\n",
    "We save the representative dataset with protein sequence and fold classification as a Pandas dataframe for further analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_json(\"../intermediate_data/foldClassification.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next step\n",
    "After you saved the dataset here, run the next step in the workflow [2-CalculateFeatures.ipynb](./2-CalculateFeatures.ipynb) or go back to [0-Workflow.ipynb](./0-Workflow.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Authors:** [Peter W. Rose](mailto:pwrose.ucsd@gmail.com), Shih-Cheng Huang, UC San Diego, October 1, 2018\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}