create_feature_matrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Feature Matrix\n",
    "This notebook serves as a guide to compile a complete feature matrix using various computational tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](figures/Fig2.png \"Title\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "import subprocess\n",
    "import itertools\n",
    "from tqdm import tqdm\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# ROOT_DIR - root directory\n",
    "ROOT_DIR = os.getcwd()+'/'\n",
    "\n",
    "# FEATURE_DIR - directory where feature dataframes are saved\n",
    "DATA_DIR = ROOT_DIR + 'dataframes/'\n",
    "\n",
    "# FASTA_DIR - directory where fasta files are saved\n",
    "FASTA_DIR = ROOT_DIR + 'fasta_files/'\n",
    "\n",
    "# MISC_DIR - directory where all other files are saved\n",
    "MISC_DIR = ROOT_DIR + 'misc_files/'\n",
    "\n",
    "# SOFTWARE_DIR - directory where all software is stored\n",
    "SOFTWARE_DIR = '/home/anand/Documents/HT_Expression/prediction_software/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of Contents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Initialize Data](#Initialize-Data)\n",
    "2. [Generate Features](#Generate-Features)\n",
    "    1. [mRNA Properties](#mRNA-Properties)\n",
    "        1. [Codon and Nucleotide Composition](#Codon-and-Nucleotide-Composition)\n",
    "        2. [tRNA Adaptation Index](#tRNA-Adaptation-Index)\n",
    "        3. [Folding Energy](#Folding-Energy)\n",
    "    2. [Physical Properties](#Physical-Properties)\n",
    "    3. [Secondary Structure](#Secondary-Structure)\n",
    "    4. [Protein Disorder](#Protein-Disorder)\n",
    "        1. [DisEMBL](#DisEMBL)\n",
    "        2. [RONN](#RONN)\n",
    "        3. [DISOPRED3](#DISOPRED3)\n",
    "3. [Compile Features](#Compile-Features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initialize Data \n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the raw dataframe into a [pandas](http://pandas.pydata.org/) dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of entries: 46521\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139  \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144  \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136  \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123  \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prest_file = DATA_DIR+'prEST_conc_data.csv'\n",
    "DF_raw_prest = pd.read_csv(prest_file,index_col=0)\n",
    "print 'Number of entries:', len(DF_raw_prest)\n",
    "DF_raw_prest.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some prESTs were expressed in duplicate experiments. Merge these results so each sequence has a single concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of total entries: 46521\n",
      "Number of unique prESTs: 45206\n"
     ]
    }
   ],
   "source": [
    "print 'Number of total entries:',len(DF_raw_prest)\n",
    "print 'Number of unique prESTs:',len(DF_raw_prest.aa_seq.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dropped = []\n",
    "DF_prest = DF_raw_prest.copy()\n",
    "for name,group in DF_raw_prest.groupby('aa_seq'):\n",
    "    if len(group)>1:\n",
    "        DF_prest.loc[group.index[0],'conc_cf'] = group.conc_cf.mean()\n",
    "        dropped += group.index[1:].tolist()\n",
    "DF_prest = DF_prest.drop(dropped).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final number of prESTs: 45206\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139  \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144  \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136  \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123  \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print 'Final number of prESTs:',len(DF_prest)\n",
    "DF_prest.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_prest.to_csv(DATA_DIR+'DF_prest.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate Features\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## mRNA Properties\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DF_mRNA_features = pd.read_csv(DATA_DIR+'DF_prest.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Directory where MFold outputs are stored\n",
    "MFOLD_DIR = MISC_DIR+'mfold_out/'\n",
    "\n",
    "# Directory that contains the tAI software\n",
    "TAI_DIR = '/home/anand/Documents/HT_Expression/prediction_software/tAI/codonR/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "# Export fasta files for all mRNA sequences and their first 40 nts\n",
    "\n",
    "sequences = []\n",
    "sequences40 = []\n",
    "\n",
    "for i,row in tqdm(DF_mRNA_features.iterrows()):\n",
    "\n",
    "    name = FASTA_DIR+'nt_sequences/'+str(row.prest_id)+'.fasta'\n",
    "    name40 = FASTA_DIR+'nt_sequences/'+str(row.prest_id)+'_40nt.fasta'\n",
    "    \n",
    "    # Export FASTA files for all mRNA sequences\n",
    "    if not os.path.isfile(name):\n",
    "        seq = SeqRecord(Seq(row.nt_seq).transcribe(), \n",
    "                        id=str(row.prest_id),description='prEST #'+str(row.prest_id))\n",
    "        sequences.append(seq)\n",
    "        with open(name,'w') as f:\n",
    "            SeqIO.write(seq,f,'fasta')\n",
    "        \n",
    "    # Export FASTA files for first 40 nts of all mRNA sequences\n",
    "    if not os.path.isfile(name40):\n",
    "        seq40 = SeqRecord(Seq(row.nt_seq).transcribe()[:40], \n",
    "                          id=str(row.prest_id),description='prEST #'+str(row.prest_id))\n",
    "        sequences40.append(seq40)\n",
    "        with open(name40,'w') as f:\n",
    "            SeqIO.write(seq40,f,'fasta')\n",
    "\n",
    "if not os.path.isfile(FASTA_DIR+'all_nt.fasta'):\n",
    "    with open(FASTA_DIR+'all_nt.fasta','w') as f:\n",
    "        SeqIO.write(sequences,f,'fasta')\n",
    "if not os.path.isfile(FASTA_DIR+'all_nt_40.fasta'):\n",
    "    with open(FASTA_DIR+'all_nt_40.fasta','w') as f:\n",
    "        SeqIO.write(sequences40,f,'fasta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Codon and Nucleotide Composition\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RNA sequences have some extra nucleotides at their head/tail, which need to be removed before selecting the codons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "0it [00:00, ?it/s]/home/anand/.virtualenvs/main/local/lib/python2.7/site-packages/Bio/Seq.py:2041: BiopythonWarning: Partial codon, len(sequence) not a multiple of three. Explicitly trim the sequence or add trailing N before translation. This may become an error in future.\n",
      "  BiopythonWarning)\n"
     ]
    }
   ],
   "source": [
    "seqs = []\n",
    "for i,row in tqdm(DF_mRNA_features.iterrows()):\n",
    "    trans_seq = Seq(row.nt_seq).translate()\n",
    "    pos = trans_seq.find(row.aa_seq)\n",
    "    seqs.append(row.nt_seq[pos*3:(pos+len(row.aa_seq))*3])\n",
    "    \n",
    "DF_mRNA_features.loc[:,'true_nt_seq'] = seqs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "# Generate a list of all codons\n",
    "\n",
    "all_codons = [''.join(a) for a in itertools.product('ATCG',repeat=3)]\n",
    "\n",
    "# For each sequence, split the sequence into groups of 3\n",
    "\n",
    "seq_codons=[]\n",
    "for i,row in tqdm(DF_mRNA_features.iterrows()):\n",
    "    seq_codons.append([row.true_nt_seq[j:j+3] for j in range(0, len(row.true_nt_seq), 3)])\n",
    "\n",
    "# Calculate the fraction of each codon appearing in each sequence\n",
    "    \n",
    "counts = [[np.true_divide(seq_codon.count(codon),len(seq_codon)) for codon in all_codons] for seq_codon in seq_codons]\n",
    "\n",
    "# Initialize the dataframe columns\n",
    "\n",
    "for codon in all_codons:\n",
    "    DF_mRNA_features[codon] = None\n",
    "    \n",
    "# Add the data to the feature matrix\n",
    "    \n",
    "DF_mRNA_features.loc[:,all_codons] = counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate nucleotide composition and GC content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for nuc in 'ATCG':\n",
    "    DF_mRNA_features.loc[:,'list_nuc_'+nuc] = [np.true_divide(seq.count(nuc),len(seq)) \\\n",
    "                                               for seq in DF_mRNA_features.true_nt_seq]\n",
    "\n",
    "DF_mRNA_features.loc[:,'GC_content'] = [np.true_divide(seq.count('G') + seq.count('C'),len(seq)) \\\n",
    "                                        for seq in DF_mRNA_features.true_nt_seq]\n",
    "DF_mRNA_features.loc[:,'GC30'] = [np.true_divide(seq[:30].count('G') + seq[:30].count('C'),len(seq[:30])) \\\n",
    "                                  for seq in DF_mRNA_features.true_nt_seq]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "      <th>true_nt_seq</th>\n",
       "      <th>AAA</th>\n",
       "      <th>AAT</th>\n",
       "      <th>AAC</th>\n",
       "      <th>...</th>\n",
       "      <th>GGA</th>\n",
       "      <th>GGT</th>\n",
       "      <th>GGC</th>\n",
       "      <th>GGG</th>\n",
       "      <th>list_nuc_A</th>\n",
       "      <th>list_nuc_T</th>\n",
       "      <th>list_nuc_C</th>\n",
       "      <th>list_nuc_G</th>\n",
       "      <th>GC_content</th>\n",
       "      <th>GC30</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "      <td>ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...</td>\n",
       "      <td>0.057554</td>\n",
       "      <td>0.043165</td>\n",
       "      <td>0.021583</td>\n",
       "      <td>...</td>\n",
       "      <td>0.007194</td>\n",
       "      <td>0.007194</td>\n",
       "      <td>0.014388</td>\n",
       "      <td>0.014388</td>\n",
       "      <td>0.306954</td>\n",
       "      <td>0.287770</td>\n",
       "      <td>0.182254</td>\n",
       "      <td>0.223022</td>\n",
       "      <td>0.405276</td>\n",
       "      <td>0.466667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "      <td>ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>...</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.027778</td>\n",
       "      <td>0.034722</td>\n",
       "      <td>0.201389</td>\n",
       "      <td>0.192130</td>\n",
       "      <td>0.275463</td>\n",
       "      <td>0.331019</td>\n",
       "      <td>0.606481</td>\n",
       "      <td>0.533333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "      <td>TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...</td>\n",
       "      <td>0.007353</td>\n",
       "      <td>0.051471</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>...</td>\n",
       "      <td>0.007353</td>\n",
       "      <td>0.014706</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.274510</td>\n",
       "      <td>0.259804</td>\n",
       "      <td>0.262255</td>\n",
       "      <td>0.203431</td>\n",
       "      <td>0.465686</td>\n",
       "      <td>0.566667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "      <td>GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>0.056911</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>...</td>\n",
       "      <td>0.024390</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.016260</td>\n",
       "      <td>0.008130</td>\n",
       "      <td>0.344173</td>\n",
       "      <td>0.205962</td>\n",
       "      <td>0.186992</td>\n",
       "      <td>0.262873</td>\n",
       "      <td>0.449864</td>\n",
       "      <td>0.366667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "      <td>CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...</td>\n",
       "      <td>0.032258</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>...</td>\n",
       "      <td>0.024194</td>\n",
       "      <td>0.016129</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.024194</td>\n",
       "      <td>0.309140</td>\n",
       "      <td>0.190860</td>\n",
       "      <td>0.247312</td>\n",
       "      <td>0.252688</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.633333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 77 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \\\n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139   \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144   \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136   \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123   \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124   \n",
       "\n",
       "                                         true_nt_seq       AAA       AAT  \\\n",
       "0  ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...  0.057554  0.043165   \n",
       "1  ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...  0.006944  0.013889   \n",
       "2  TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...  0.007353  0.051471   \n",
       "3  GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...  0.040650  0.056911   \n",
       "4  CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...  0.032258  0.008065   \n",
       "\n",
       "        AAC    ...          GGA       GGT       GGC       GGG  list_nuc_A  \\\n",
       "0  0.021583    ...     0.007194  0.007194  0.014388  0.014388    0.306954   \n",
       "1  0.013889    ...     0.006944  0.006944  0.027778  0.034722    0.201389   \n",
       "2  0.058824    ...     0.007353  0.014706  0.000000  0.000000    0.274510   \n",
       "3  0.040650    ...     0.024390  0.000000  0.016260  0.008130    0.344173   \n",
       "4  0.008065    ...     0.024194  0.016129  0.000000  0.024194    0.309140   \n",
       "\n",
       "   list_nuc_T  list_nuc_C  list_nuc_G  GC_content      GC30  \n",
       "0    0.287770    0.182254    0.223022    0.405276  0.466667  \n",
       "1    0.192130    0.275463    0.331019    0.606481  0.533333  \n",
       "2    0.259804    0.262255    0.203431    0.465686  0.566667  \n",
       "3    0.205962    0.186992    0.262873    0.449864  0.366667  \n",
       "4    0.190860    0.247312    0.252688    0.500000  0.633333  \n",
       "\n",
       "[5 rows x 77 columns]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DF_mRNA_features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SD Sequences\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count the number of Shine-Delgarno (SD) and SD-like sequences on both the forward and reverse strands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def complement(seq):\n",
    "    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} \n",
    "    bases = list(seq) \n",
    "    bases = [complement[base] for base in bases] \n",
    "    return ''.join(bases)\n",
    "\n",
    "def reverse_complement(s):\n",
    "    return complement(s[::-1])\n",
    "\n",
    "def count_sd_like(seq):\n",
    "    sd_like_pat=r'(?=(AGG)|(GGA)|(GAG)|(GGG)|(GGT)|(GTG))'\n",
    "    revcomp = reverse_complement(seq)\n",
    "    return [len([m.start() for m in re.finditer(sd_like_pat, seq)]),\n",
    "            len([m.start() for m in re.finditer(sd_like_pat, revcomp)])]\n",
    "\n",
    "def count_sd(seq):\n",
    "    sd_pat = r'(?=(GG.GG.))'\n",
    "    revcomp = reverse_complement(seq)\n",
    "    return [len([m.start() for m in re.finditer(sd_pat, seq)]),\n",
    "            len([m.start() for m in re.finditer(sd_pat, revcomp)])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sd_like = [count_sd_like(seq) for seq in DF_mRNA_features.true_nt_seq]\n",
    "\n",
    "DF_mRNA_features['sd_like_fwd'] = [num[0] for num in sd_like]\n",
    "DF_mRNA_features['sd_like_rev'] = [num[1] for num in sd_like]\n",
    "DF_mRNA_features['sd_like_fwd_frac'] = np.true_divide(DF_mRNA_features['sd_like_fwd'],\n",
    "                                                      [len(seq) for seq in DF_mRNA_features.true_nt_seq])\n",
    "DF_mRNA_features['sd_like_rev_frac'] = np.true_divide(DF_mRNA_features['sd_like_rev'],\n",
    "                                                      [len(seq) for seq in DF_mRNA_features.true_nt_seq])\n",
    "\n",
    "\n",
    "sd = [count_sd(seq) for seq in DF_mRNA_features.true_nt_seq]\n",
    "DF_mRNA_features['sd_seq_fwd'] = [num[0] for num in sd]\n",
    "DF_mRNA_features['sd_seq_rev'] = [num[1] for num in sd]\n",
    "DF_mRNA_features['sd_seq_fwd_frac'] = np.true_divide(DF_mRNA_features['sd_seq_fwd'],\n",
    "                                                     [len(seq) for seq in DF_mRNA_features.true_nt_seq])\n",
    "DF_mRNA_features['sd_seq_rev_frac'] = np.true_divide(DF_mRNA_features['sd_seq_rev'],\n",
    "                                                     [len(seq) for seq in DF_mRNA_features.true_nt_seq])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "      <th>true_nt_seq</th>\n",
       "      <th>AAA</th>\n",
       "      <th>AAT</th>\n",
       "      <th>AAC</th>\n",
       "      <th>...</th>\n",
       "      <th>GC_content</th>\n",
       "      <th>GC30</th>\n",
       "      <th>sd_like_fwd</th>\n",
       "      <th>sd_like_rev</th>\n",
       "      <th>sd_like_fwd_frac</th>\n",
       "      <th>sd_like_rev_frac</th>\n",
       "      <th>sd_seq_fwd</th>\n",
       "      <th>sd_seq_rev</th>\n",
       "      <th>sd_seq_fwd_frac</th>\n",
       "      <th>sd_seq_rev_frac</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "      <td>ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...</td>\n",
       "      <td>0.057554</td>\n",
       "      <td>0.043165</td>\n",
       "      <td>0.021583</td>\n",
       "      <td>...</td>\n",
       "      <td>0.405276</td>\n",
       "      <td>0.466667</td>\n",
       "      <td>32</td>\n",
       "      <td>28</td>\n",
       "      <td>0.076739</td>\n",
       "      <td>0.067146</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0.004796</td>\n",
       "      <td>0.002398</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "      <td>ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>...</td>\n",
       "      <td>0.606481</td>\n",
       "      <td>0.533333</td>\n",
       "      <td>77</td>\n",
       "      <td>56</td>\n",
       "      <td>0.178241</td>\n",
       "      <td>0.129630</td>\n",
       "      <td>10</td>\n",
       "      <td>3</td>\n",
       "      <td>0.023148</td>\n",
       "      <td>0.006944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "      <td>TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...</td>\n",
       "      <td>0.007353</td>\n",
       "      <td>0.051471</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>...</td>\n",
       "      <td>0.465686</td>\n",
       "      <td>0.566667</td>\n",
       "      <td>34</td>\n",
       "      <td>52</td>\n",
       "      <td>0.083333</td>\n",
       "      <td>0.127451</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>0.002451</td>\n",
       "      <td>0.009804</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "      <td>GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>0.056911</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>...</td>\n",
       "      <td>0.449864</td>\n",
       "      <td>0.366667</td>\n",
       "      <td>40</td>\n",
       "      <td>21</td>\n",
       "      <td>0.108401</td>\n",
       "      <td>0.056911</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0.005420</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "      <td>CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...</td>\n",
       "      <td>0.032258</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>...</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.633333</td>\n",
       "      <td>41</td>\n",
       "      <td>29</td>\n",
       "      <td>0.110215</td>\n",
       "      <td>0.077957</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.002688</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 85 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \\\n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139   \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144   \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136   \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123   \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124   \n",
       "\n",
       "                                         true_nt_seq       AAA       AAT  \\\n",
       "0  ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...  0.057554  0.043165   \n",
       "1  ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...  0.006944  0.013889   \n",
       "2  TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...  0.007353  0.051471   \n",
       "3  GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...  0.040650  0.056911   \n",
       "4  CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...  0.032258  0.008065   \n",
       "\n",
       "        AAC       ...         GC_content      GC30  sd_like_fwd  sd_like_rev  \\\n",
       "0  0.021583       ...           0.405276  0.466667           32           28   \n",
       "1  0.013889       ...           0.606481  0.533333           77           56   \n",
       "2  0.058824       ...           0.465686  0.566667           34           52   \n",
       "3  0.040650       ...           0.449864  0.366667           40           21   \n",
       "4  0.008065       ...           0.500000  0.633333           41           29   \n",
       "\n",
       "   sd_like_fwd_frac  sd_like_rev_frac  sd_seq_fwd  sd_seq_rev  \\\n",
       "0          0.076739          0.067146           2           1   \n",
       "1          0.178241          0.129630          10           3   \n",
       "2          0.083333          0.127451           1           4   \n",
       "3          0.108401          0.056911           2           0   \n",
       "4          0.110215          0.077957           1           0   \n",
       "\n",
       "   sd_seq_fwd_frac  sd_seq_rev_frac  \n",
       "0         0.004796         0.002398  \n",
       "1         0.023148         0.006944  \n",
       "2         0.002451         0.009804  \n",
       "3         0.005420         0.000000  \n",
       "4         0.002688         0.000000  \n",
       "\n",
       "[5 rows x 85 columns]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DF_mRNA_features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### tRNA Adaptation Index\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the tAI code on all 45206 prESTs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "command1 = ['perl', TAI_DIR+'codonM', FASTA_DIR+'all_nt.fasta',MISC_DIR+'tAI_files/prEST_tAI.m']\n",
    "print 'Running Command:'+' '.join(command1) + '\\n'\n",
    "\n",
    "print subprocess.check_output(command1)\n",
    "\n",
    "command2 = ['Rscript',MISC_DIR+'tAI_files/calc_tAI.R',MISC_DIR+'tAI_files/',TAI_DIR]\n",
    "print 'Running Command:'+' '.join(command2) + '\\n'\n",
    "\n",
    "print subprocess.check_output(command2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now add the tAI as a feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DF_mRNA_features.loc[:,'tAI'] = pd.read_csv(MISC_DIR+'tAI_files/prEST_tAI.csv')['x']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Folding Energy\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculation of all folding energies may take up to 15 hours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "  0%|          | 0/45206 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Command: 'mfold SEQ='/home/anand/Documents/HT_Expression/ML_HPA/fasta_files/nt_sequences/PREST_ID.fasta' MAX=1'\n"
     ]
    }
   ],
   "source": [
    "print 'Command: \\'mfold SEQ=\\''+FASTA_DIR+'nt_sequences/PREST_ID.fasta\\' MAX=1\\''\n",
    "\n",
    "os.chdir(MFOLD_DIR)\n",
    "\n",
    "energy = []\n",
    "for prest_id in tqdm(DF_mRNA_features.prest_id):\n",
    "\n",
    "    # Run mRNA folding energy prediction software on full mRNA\n",
    "    if not os.path.isfile(MFOLD_DIR+str(prest_id)+'.ct'):\n",
    "        command = 'mfold SEQ=\\''+FASTA_DIR+'nt_sequences/'+str(prest_id)+'.fasta\\' MAX=1'\n",
    "        os.system(command)\n",
    "    \n",
    "    # Run mRNA folding energy prediction software on first 40 nucleotides\n",
    "    if not os.path.isfile(MFOLD_DIR+str(prest_id)+'_40nt'+'.ct'):\n",
    "        command = 'mfold SEQ=\\''+FASTA_DIR+'nt_sequences/'+str(prest_id)+'_40nt.fasta\\' MAX=1'\n",
    "        os.system(command)\n",
    "    \n",
    "    # Remove non-essential files\n",
    "    for f in os.listdir(MFOLD_DIR):\n",
    "        if not f.endswith('.ct') or f.endswith('_1.ct'):\n",
    "            os.remove(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse the output of the folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "for i,row in tqdm(DF_mRNA_features.iterrows()):\n",
    "\n",
    "    if os.path.isfile(MFOLD_DIR+str(row.prest_id)+'.ct'):\n",
    "        with open(MFOLD_DIR+str(row.prest_id)+'.ct','r') as f:\n",
    "            result = f.read()\n",
    "            energy = re.search('-\\d*\\.\\d*',result)\n",
    "            if not energy:\n",
    "                energy = re.search('\\d*\\.\\d*',result)\n",
    "            DF_mRNA_features.loc[i,'RNA_folding_energy'] = float(energy.group())\n",
    "            \n",
    "    if os.path.isfile(MFOLD_DIR+str(row.prest_id)+'_40nt.ct'):\n",
    "        with open(MFOLD_DIR+str(row.prest_id)+'_40nt.ct','r') as f:\n",
    "            result = f.read()\n",
    "            energy = re.search('-\\d*\\.\\d*',result)\n",
    "            if not energy:\n",
    "                energy = re.search('\\d*\\.\\d*',result)\n",
    "            DF_mRNA_features.loc[i,'RNA_40_energy'] = float(energy.group())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "      <th>true_nt_seq</th>\n",
       "      <th>AAA</th>\n",
       "      <th>AAT</th>\n",
       "      <th>AAC</th>\n",
       "      <th>...</th>\n",
       "      <th>sd_like_rev</th>\n",
       "      <th>sd_like_fwd_frac</th>\n",
       "      <th>sd_like_rev_frac</th>\n",
       "      <th>sd_seq_fwd</th>\n",
       "      <th>sd_seq_rev</th>\n",
       "      <th>sd_seq_fwd_frac</th>\n",
       "      <th>sd_seq_rev_frac</th>\n",
       "      <th>tAI</th>\n",
       "      <th>RNA_folding_energy</th>\n",
       "      <th>RNA_40_energy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "      <td>ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...</td>\n",
       "      <td>0.057554</td>\n",
       "      <td>0.043165</td>\n",
       "      <td>0.021583</td>\n",
       "      <td>...</td>\n",
       "      <td>28</td>\n",
       "      <td>0.076739</td>\n",
       "      <td>0.067146</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0.004796</td>\n",
       "      <td>0.002398</td>\n",
       "      <td>0.308574</td>\n",
       "      <td>-111.56</td>\n",
       "      <td>-3.70</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "      <td>ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>0.013889</td>\n",
       "      <td>...</td>\n",
       "      <td>56</td>\n",
       "      <td>0.178241</td>\n",
       "      <td>0.129630</td>\n",
       "      <td>10</td>\n",
       "      <td>3</td>\n",
       "      <td>0.023148</td>\n",
       "      <td>0.006944</td>\n",
       "      <td>0.252830</td>\n",
       "      <td>-192.52</td>\n",
       "      <td>-10.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "      <td>TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...</td>\n",
       "      <td>0.007353</td>\n",
       "      <td>0.051471</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>...</td>\n",
       "      <td>52</td>\n",
       "      <td>0.083333</td>\n",
       "      <td>0.127451</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>0.002451</td>\n",
       "      <td>0.009804</td>\n",
       "      <td>0.230369</td>\n",
       "      <td>-113.76</td>\n",
       "      <td>-4.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "      <td>GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>0.056911</td>\n",
       "      <td>0.040650</td>\n",
       "      <td>...</td>\n",
       "      <td>21</td>\n",
       "      <td>0.108401</td>\n",
       "      <td>0.056911</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0.005420</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.356271</td>\n",
       "      <td>-97.22</td>\n",
       "      <td>-5.74</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "      <td>CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...</td>\n",
       "      <td>0.032258</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>0.008065</td>\n",
       "      <td>...</td>\n",
       "      <td>29</td>\n",
       "      <td>0.110215</td>\n",
       "      <td>0.077957</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.002688</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.296122</td>\n",
       "      <td>-107.98</td>\n",
       "      <td>-9.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 88 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \\\n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139   \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144   \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136   \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123   \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124   \n",
       "\n",
       "                                         true_nt_seq       AAA       AAT  \\\n",
       "0  ATTATGACAGCTCCCTCCAGTTTTGAGCAGTTTAAAGTGGCAATGA...  0.057554  0.043165   \n",
       "1  ACCTACTATGCCTGGAAGCATGAGCTGCTGGGCTCTGGCACCTGCC...  0.006944  0.013889   \n",
       "2  TCACTCCATGCCAGACCCCCACAGTTTACGAGGGCTCAGTGGTTTG...  0.007353  0.051471   \n",
       "3  GCGAGAGCATTAAATGAAAGCAAAAGAGTTAATAATGGCAACACGG...  0.040650  0.056911   \n",
       "4  CATCGGAAAGAGCCTGGGGCAAGGCTGGAGGCCACAAGAGGAGCTG...  0.032258  0.008065   \n",
       "\n",
       "        AAC      ...        sd_like_rev  sd_like_fwd_frac  sd_like_rev_frac  \\\n",
       "0  0.021583      ...                 28          0.076739          0.067146   \n",
       "1  0.013889      ...                 56          0.178241          0.129630   \n",
       "2  0.058824      ...                 52          0.083333          0.127451   \n",
       "3  0.040650      ...                 21          0.108401          0.056911   \n",
       "4  0.008065      ...                 29          0.110215          0.077957   \n",
       "\n",
       "   sd_seq_fwd  sd_seq_rev  sd_seq_fwd_frac  sd_seq_rev_frac       tAI  \\\n",
       "0           2           1         0.004796         0.002398  0.308574   \n",
       "1          10           3         0.023148         0.006944  0.252830   \n",
       "2           1           4         0.002451         0.009804  0.230369   \n",
       "3           2           0         0.005420         0.000000  0.356271   \n",
       "4           1           0         0.002688         0.000000  0.296122   \n",
       "\n",
       "   RNA_folding_energy  RNA_40_energy  \n",
       "0             -111.56          -3.70  \n",
       "1             -192.52         -10.00  \n",
       "2             -113.76          -4.00  \n",
       "3              -97.22          -5.74  \n",
       "4             -107.98          -9.00  \n",
       "\n",
       "[5 rows x 88 columns]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DF_mRNA_features.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_mRNA_features.to_csv(DATA_DIR+'mRNA_features.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Physical Properties\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_physical_features = pd.read_csv(DATA_DIR+'DF_prest.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BioPython ProtParam quickly calculates many physical properties for any peptide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 45206/45206 [00:12<00:00, 3730.07it/s]\n"
     ]
    }
   ],
   "source": [
    "from Bio.SeqUtils import ProtParam\n",
    "\n",
    "mw = []\n",
    "pI= []\n",
    "arom = []\n",
    "instab = []\n",
    "gravy = []\n",
    "for seq in tqdm(DF_physical_features.aa_seq):\n",
    "    data = ProtParam.ProteinAnalysis(seq)\n",
    "    mw.append(data.molecular_weight())\n",
    "    pI.append(data.isoelectric_point())\n",
    "    arom.append(data.aromaticity())\n",
    "    instab.append(data.instability_index())\n",
    "    gravy.append(data.gravy())\n",
    "    \n",
    "DF_physical_features['bio_pI'] = pI\n",
    "DF_physical_features['bio_mW'] = mw\n",
    "DF_physical_features['bio_aromaticity'] = arom\n",
    "DF_physical_features['bio_instability'] = instab\n",
    "DF_physical_features['bio_gravy'] = gravy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the counts of each amino acid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "amino_acids = 'ARNDCQEGHILKMFPSTWYV'\n",
    "counts = [[np.true_divide(seq.count(aa),len(seq)) for seq in DF_physical_features.aa_seq] for aa in amino_acids]\n",
    "for i in range(len(amino_acids)):\n",
    "    DF_physical_features['list_comp_'+amino_acids[i]] = counts[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the charge of each prEST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_physical_features['charge'] = [seq.count('K')+seq.count('R')-seq.count('D')-seq.count('E') for seq in DF_physical_features.aa_seq]\n",
    "DF_physical_features['abs_charge']=abs(DF_physical_features['charge'])\n",
    "DF_physical_features['avg_charge']=np.true_divide(DF_physical_features['charge'],[len(seq) for seq in DF_physical_features['aa_seq']])\n",
    "DF_physical_features['abs_avg_charge'] = abs(DF_physical_features['avg_charge'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the presence of various amino acid types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_physical_features['frac_aliphatic'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'AGILPV'])\n",
    "DF_physical_features['frac_uncharged_polar'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'STNQ'])\n",
    "DF_physical_features['frac_polar'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'QNHSTYCMW'])\n",
    "DF_physical_features['frac_hydrophobic'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'AGILPVF'])\n",
    "DF_physical_features['frac_positive'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'HKR'])\n",
    "DF_physical_features['frac_sulfur'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'CM'])\n",
    "DF_physical_features['frac_negative'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'DE'])\n",
    "DF_physical_features['frac_amide'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'NQ'])\n",
    "DF_physical_features['frac_alcohol'] = sum([DF_physical_features['list_comp_'+aa] for aa in 'ST'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_physical_features.to_csv(DATA_DIR+'physical_features.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Secondary Structure\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SCRATCH-1D predictor is computationally intensive and may take weeks to run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_ss_features = pd.read_csv(DATA_DIR+'DF_prest.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "# Export fasta files for all sequences\n",
    "\n",
    "sequences = []\n",
    "\n",
    "for i,row in tqdm(DF_ss_features.iterrows()):\n",
    "\n",
    "    name = FASTA_DIR+'aa_sequences/'+str(row.prest_id)+'.fasta'\n",
    "    \n",
    "    # Export FASTA files for all mRNA sequences\n",
    "    if not os.path.isfile(name):\n",
    "        seq = SeqRecord(Seq(row.aa_seq), id=str(row.prest_id),description='prEST #'+str(row.prest_id))\n",
    "        sequences.append(seq)\n",
    "        with open(name,'w') as f:\n",
    "            SeqIO.write(seq,f,'fasta')\n",
    "        \n",
    "if not os.path.isfile(FASTA_DIR+'all_prEST.fasta'):\n",
    "    with open(FASTA_DIR+'all_prEST.fasta','w') as f:\n",
    "        SeqIO.write(sequences,f,'fasta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running command: /home/anand/Documents/HT_Expression/prediction_software/SCRATCH-1D_1.1/bin/run_SCRATCH-1D_predictors.sh /home/anand/Documents/HT_Expression/ML_HPA/fasta_files/all_prEST.fasta /home/anand/Documents/HT_Expression/ML_HPA/misc_files/scratch_out/all_prEST 8\n"
     ]
    }
   ],
   "source": [
    "scratch_cmd = '/home/anand/Documents/HT_Expression/prediction_software/SCRATCH-1D_1.1/bin/run_SCRATCH-1D_predictors.sh'\n",
    "name = FASTA_DIR+'all_prEST.fasta'\n",
    "out = MISC_DIR+'scratch_out/all_prEST'\n",
    "n_cores = '8'\n",
    "\n",
    "print 'Running command:',scratch_cmd,name,out,n_cores\n",
    "\n",
    "if not os.path.isfile(out+'.ss'):\n",
    "    subprocess.call([scratch_cmd,name,out,n_cores])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SCRATCH-1D predictor returns 4 outputs:\n",
    "* .ss: 3-letter secondary structure\n",
    "* .ss8: 8-letter secondary structure\n",
    "* .acc: Solvent accessible residues\n",
    "* .acc20: Solvent accessibility index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## SS\n",
    "\n",
    "ss_list = []\n",
    "for record in SeqIO.parse(MISC_DIR+'scratch_out/all_prEST.ss','fasta'):\n",
    "    ss_list.append([int(record.id),str(record.seq)])\n",
    "df_ss = pd.DataFrame(ss_list,columns=['prest_id','SCRATCH_ss'])\n",
    "DF_ss_features = pd.merge(DF_ss_features,df_ss,on='prest_id')\n",
    "\n",
    "## SS8\n",
    "ss8_list = []\n",
    "for record in SeqIO.parse(MISC_DIR+'scratch_out/all_prEST.ss8','fasta'):\n",
    "    ss8_list.append([int(record.id),str(record.seq)])\n",
    "df_ss8 = pd.DataFrame(ss8_list,columns=['prest_id','SCRATCH_ss8'])\n",
    "DF_ss_features = pd.merge(DF_ss_features,df_ss8,on='prest_id')\n",
    "\n",
    "## ACC\n",
    "acc_list = []\n",
    "for record in SeqIO.parse(MISC_DIR+'scratch_out/all_prEST.acc','fasta'):\n",
    "    acc_list.append([int(record.id),str(record.seq)])\n",
    "df_acc = pd.DataFrame(acc_list,columns=['prest_id','SCRATCH_acc'])\n",
    "DF_ss_features = pd.merge(DF_ss_features,df_acc,on='prest_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## ACC20\n",
    "acc20_list = []\n",
    "with open(MISC_DIR+'scratch_out/all_prEST.acc20','r') as f:\n",
    "    result = f.read()\n",
    "    lines = result[1:-1].split('\\n>')\n",
    "    for line in lines:\n",
    "        [name,nums] = line.split('\\n')\n",
    "        acc20_mean = np.mean([int(a) for a in nums.split(' ')])\n",
    "        seq_id = int(name[:name.find(' ')])\n",
    "        acc20_list.append([seq_id,acc20_mean])\n",
    "df_acc20 = pd.DataFrame(acc20_list,columns=['prest_id','acc20_mean'])\n",
    "DF_ss_features = pd.merge(DF_ss_features,df_acc20,on='prest_id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct features based on the SCRATCH-1D output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ss_names = ['ss_helix','ss_ext','ss_c']\n",
    "ss_letters = 'HEC'\n",
    "ss8_names = ['ss8_helix','ss8_ext','ss8_turn','ss8_helix3',\n",
    "             'ss8_pi_helix','ss8_bridge','ss8_bend','ss8_coil']\n",
    "ss8_letters = 'HETGIBSC'\n",
    "\n",
    "for i in range(3):\n",
    "    DF_ss_features[ss_names[i]] = [np.true_divide(ss.count(ss_letters[i]),len(ss)) \\\n",
    "                                   for ss in DF_ss_features.SCRATCH_ss]\n",
    "\n",
    "for i in range(8):\n",
    "    DF_ss_features[ss8_names[i]] = [np.true_divide(ss.count(ss8_letters[i]),len(ss)) \\\n",
    "                                    for ss in DF_ss_features.SCRATCH_ss8]\n",
    "    \n",
    "DF_ss_features['acc_frac'] = [np.true_divide(seq.count('e'),len(seq)) \\\n",
    "                              for seq in DF_ss_features.SCRATCH_acc]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate fraction of hydrophilic and hydrophobic solvent accessible and inaccessible residues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "hydrophilic='RDEKSNQ'\n",
    "hydrophobic='WFYLIVMACHT'\n",
    "\n",
    "for i,row in tqdm(DF_ss_features.iterrows()):\n",
    "    acc_hydrophobic_out = 0\n",
    "    acc_hydrophilic_out = 0\n",
    "    acc_hydrophobic_in = 0\n",
    "    acc_hydrophilic_in = 0\n",
    "    for j in xrange(len(row.aa_seq)):\n",
    "        if row.aa_seq[j] in hydrophobic:\n",
    "            if row.SCRATCH_acc[j]=='e':\n",
    "                acc_hydrophobic_out+=1\n",
    "            else:\n",
    "                acc_hydrophobic_in+=1\n",
    "        elif row.aa_seq[j] in hydrophilic:\n",
    "            if row.SCRATCH_acc[j]=='e':\n",
    "                acc_hydrophilic_out+=1\n",
    "            else:\n",
    "                acc_hydrophilic_in+=1\n",
    "    DF_ss_features.loc[i,'acc_hydrophobic_out'] = np.true_divide(acc_hydrophobic_out,row.aa_len)\n",
    "    DF_ss_features.loc[i,'acc_hydrophobic_in'] = np.true_divide(acc_hydrophobic_in,row.aa_len)\n",
    "    DF_ss_features.loc[i,'acc_hydrophilic_out'] = np.true_divide(acc_hydrophilic_out,row.aa_len)\n",
    "    DF_ss_features.loc[i,'acc_hydrophilic_in'] = np.true_divide(acc_hydrophilic_in,row.aa_len)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the average GRAVY of solvent accessible and inaccessible residues.\n",
    "\n",
    "kd is a dictionary that contains the Kyle-Doolittle index of hydrophobicity used to calculate GRAVY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "from Bio.SeqUtils.ProtParamData import kd\n",
    "\n",
    "for i,row in tqdm(DF_ss_features.iterrows()):\n",
    "    out_gravy = 0\n",
    "    in_gravy = 0\n",
    "    for j in xrange(len(row.aa_seq)):\n",
    "        if row.SCRATCH_acc[j] == 'e':\n",
    "            out_gravy += kd[row.aa_seq[j]]\n",
    "        else:\n",
    "            in_gravy += kd[row.aa_seq[j]]\n",
    "    DF_ss_features.loc[i,'out_gravy'] = np.true_divide(out_gravy,len(row.aa_seq))\n",
    "    DF_ss_features.loc[i,'in_gravy'] = np.true_divide(in_gravy,len(row.aa_seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_ss_features.to_csv(DATA_DIR+'secondary_structure.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Protein Disorder\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_disorder_features = pd.read_csv(DATA_DIR+'DF_prest.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DisEMBL\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run DisEMBL on all prESTs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "disembl_cmd = '/home/anand/Documents/HT_Expression/prediction_software/DisEMBL-1.4/DisEMBL.py'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "DF_disorder_features['disembl_COILS'] = None\n",
    "DF_disorder_features['disembl_REM465'] = None\n",
    "DF_disorder_features['disembl_HOTLOOPS'] = None\n",
    "\n",
    "disembl_list = []\n",
    "\n",
    "for prest_id in tqdm(DF_disorder_features.prest_id):\n",
    "    out = subprocess.check_output([disembl_cmd,'8','8','4','1.2','1.4','1.2',\n",
    "                                   FASTA_DIR+'aa_sequences/'+str(prest_id)+'.fasta']).split('\\n')\n",
    "    disembl_list.append([out[9],out[11],out[13]])\n",
    "    \n",
    "DF_disorder_features[['disembl_COILS','disembl_REM465','disembl_HOTLOOPS']] = disembl_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import string\n",
    "\n",
    "def count_upper(seq):\n",
    "    return len(filter(lambda x: x in string.uppercase, seq))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DF_disorder_features.loc[:,'disembl_COILS_frac'] = [np.true_divide(count_upper(seq),len(seq)) for seq in DF_disorder_features.disembl_COILS]\n",
    "DF_disorder_features.loc[:,'disembl_REM465_frac'] = [np.true_divide(count_upper(seq),len(seq)) for seq in DF_disorder_features.disembl_REM465]\n",
    "DF_disorder_features.loc[:,'disembl_HOTLOOPS_frac'] = [np.true_divide(count_upper(seq),len(seq)) for seq in DF_disorder_features.disembl_HOTLOOPS]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RONN\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RONN_DIR = '/home/anand/Documents/HT_Expression/prediction_software/RONN'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run RONN Command on all prESTs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if not os.path.isfile(MISC_DIR+'RONN_out/out.txt'):\n",
    "    print 'Running command:',' '.join(['bash',MISC_DIR+'RONN_out/ronn_script.sh',FASTA_DIR+'aa_sequences',RONN_DIR])\n",
    "    subprocess.call(['bash',MISC_DIR+'RONN_out/ronn_script.sh',FASTA_DIR+'aa_sequences',RONN_DIR])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "ronn_results = {}\n",
    "ronn_avg={}\n",
    "with open(MISC_DIR+'RONN_out/out.txt', 'rb') as f:\n",
    "    entries = f.read().split('>')[1:]\n",
    "    for seq in tqdm(entries):\n",
    "        seq_list = seq.split('\\n')\n",
    "        title = seq_list[0][:seq_list[0].find(' ')]\n",
    "        disorder = [aa.split('\\t') for aa in seq_list[1:] if aa!='']\n",
    "        disorder = [[aa[0],float(aa[1])] for aa in disorder]\n",
    "        ronn_avg[title] = np.mean([aa[1] for aa in disorder])\n",
    "        ronn_results[title] = ''.join([aa[0].upper() if aa[1]>.5 else aa[0].lower() for aa in disorder])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the features to the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "for i,row in tqdm(DF_disorder_features.iterrows()):\n",
    "    DF_disorder_features.loc[i,'ronn_avg'] = ronn_avg[str(row.prest_id)]\n",
    "    DF_disorder_features.loc[i,'ronn_results'] = ronn_results[str(row.prest_id)]\n",
    "\n",
    "DF_disorder_features['ronn_frac'] = [np.true_divide(count_upper(seq),len(seq)) for seq in DF_disorder_features.ronn_results]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prest_id</th>\n",
       "      <th>uniprot_id</th>\n",
       "      <th>conc_cf</th>\n",
       "      <th>aa_seq</th>\n",
       "      <th>nt_seq</th>\n",
       "      <th>aa_len</th>\n",
       "      <th>disembl_COILS</th>\n",
       "      <th>disembl_REM465</th>\n",
       "      <th>disembl_HOTLOOPS</th>\n",
       "      <th>disembl_COILS_frac</th>\n",
       "      <th>...</th>\n",
       "      <th>disembl_HOTLOOPS_frac</th>\n",
       "      <th>ronn_avg</th>\n",
       "      <th>ronn_results</th>\n",
       "      <th>ronn_frac</th>\n",
       "      <th>disopred_results</th>\n",
       "      <th>disopred_avg</th>\n",
       "      <th>disopred_pb_results</th>\n",
       "      <th>disopred_pb_avg</th>\n",
       "      <th>disopred_frac</th>\n",
       "      <th>disopred_pb_frac</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>140095</td>\n",
       "      <td>G3V3N0</td>\n",
       "      <td>4.3075</td>\n",
       "      <td>IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...</td>\n",
       "      <td>139</td>\n",
       "      <td>imtapssfeqfkvamnylqLYNVPDCLEDIQDADCSSSKCSSSASS...</td>\n",
       "      <td>imtapssfeqfkvamnylqlynvpdclediqDADCSSSKCSSSASS...</td>\n",
       "      <td>IMTAPSSFeqfkvamnylqlynvpdclediqdadcSSSKCSSSASS...</td>\n",
       "      <td>0.769784</td>\n",
       "      <td>...</td>\n",
       "      <td>0.410072</td>\n",
       "      <td>0.495827</td>\n",
       "      <td>imtapssfeqfkvamnylqlynvpdclediqDADCSSSKCSSSASS...</td>\n",
       "      <td>0.482014</td>\n",
       "      <td>ImtapssfeqfkvamnylqlynvpdclediqdadcsSSKCSSSASS...</td>\n",
       "      <td>0.387050</td>\n",
       "      <td>-...................................^^^^^^^^^^...</td>\n",
       "      <td>0.317770</td>\n",
       "      <td>0.388489</td>\n",
       "      <td>0.359712</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>140099</td>\n",
       "      <td>G3V537</td>\n",
       "      <td>2.9154</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...</td>\n",
       "      <td>144</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQvaeeelecsa...</td>\n",
       "      <td>tyyawkhellgsgtcpalpprevlgmeeleklpEEQVAEEELecsa...</td>\n",
       "      <td>TYYAWKHELLGSGTCPALPPrevlgmeeleklpeeqvaeeelecsa...</td>\n",
       "      <td>0.791667</td>\n",
       "      <td>...</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.488125</td>\n",
       "      <td>tyyawkhellGSGTcPALPPREVLGMEELEKLPEEQVAEEELECSa...</td>\n",
       "      <td>0.493056</td>\n",
       "      <td>tyyawkhellgsgtCPALPPREVLGMEELEKLPEEQVAEEELECSA...</td>\n",
       "      <td>0.421597</td>\n",
       "      <td>..............^^^^^^--^^^^^^^^^^^^^^^^^^^-^^^^...</td>\n",
       "      <td>0.352917</td>\n",
       "      <td>0.465278</td>\n",
       "      <td>0.423611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>140225</td>\n",
       "      <td>P12724</td>\n",
       "      <td>1.4877</td>\n",
       "      <td>SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...</td>\n",
       "      <td>136</td>\n",
       "      <td>SLHARPPQFTraqwfaiqHISLNPPRCTiamrainnyrwrcknqnt...</td>\n",
       "      <td>slharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...</td>\n",
       "      <td>SLHARPPQftraqwfaiqhislnpprctiamrainnyrwrcknqnt...</td>\n",
       "      <td>0.720588</td>\n",
       "      <td>...</td>\n",
       "      <td>0.338235</td>\n",
       "      <td>0.368750</td>\n",
       "      <td>SLharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...</td>\n",
       "      <td>0.058824</td>\n",
       "      <td>slharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...</td>\n",
       "      <td>0.034338</td>\n",
       "      <td>.................................................</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>140235</td>\n",
       "      <td>H0YH02</td>\n",
       "      <td>6.7224</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...</td>\n",
       "      <td>123</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>0.617886</td>\n",
       "      <td>...</td>\n",
       "      <td>0.552846</td>\n",
       "      <td>0.571382</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>0.536585</td>\n",
       "      <td>ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...</td>\n",
       "      <td>0.425366</td>\n",
       "      <td>^^^^^^^^^^^^^^^-^---^-^^^^^^^^^^^^^--^^^^^^^^^...</td>\n",
       "      <td>0.295122</td>\n",
       "      <td>0.430894</td>\n",
       "      <td>0.373984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>140309</td>\n",
       "      <td>F5GYC5</td>\n",
       "      <td>3.3848</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...</td>\n",
       "      <td>124</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMitrpSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQlgegkcpssqhl...</td>\n",
       "      <td>0.693548</td>\n",
       "      <td>...</td>\n",
       "      <td>0.274194</td>\n",
       "      <td>0.556613</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...</td>\n",
       "      <td>0.620968</td>\n",
       "      <td>HRKEPGARLEATRGAARPHKQGTkpmitrpsvsqlgegkcpssqhl...</td>\n",
       "      <td>0.236613</td>\n",
       "      <td>^^^^^^^^^^^^^^^^^^^^^^^^.........................</td>\n",
       "      <td>0.185726</td>\n",
       "      <td>0.201613</td>\n",
       "      <td>0.209677</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 21 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   prest_id uniprot_id  conc_cf  \\\n",
       "0    140095     G3V3N0   4.3075   \n",
       "1    140099     G3V537   2.9154   \n",
       "2    140225     P12724   1.4877   \n",
       "3    140235     H0YH02   6.7224   \n",
       "4    140309     F5GYC5   3.3848   \n",
       "\n",
       "                                              aa_seq  \\\n",
       "0  IMTAPSSFEQFKVAMNYLQLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQVAEEELECSA...   \n",
       "2  SLHARPPQFTRAQWFAIQHISLNPPRCTIAMRAINNYRWRCKNQNT...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                              nt_seq  aa_len  \\\n",
       "0  GACAAGCTTGCGGCCGCAATTATGACAGCTCCCTCCAGTTTTGAGC...     139   \n",
       "1  GACAAGCTTGCGGCCGCAACCTACTATGCCTGGAAGCATGAGCTGC...     144   \n",
       "2  GACAAGCTTGCGGCCGCATCACTCCATGCCAGACCCCCACAGTTTA...     136   \n",
       "3  GACAAGCTTGCGGCCGCAGCGAGAGCATTAAATGAAAGCAAAAGAG...     123   \n",
       "4  GACAAGCTTGCGGCCGCACATCGGAAAGAGCCTGGGGCAAGGCTGG...     124   \n",
       "\n",
       "                                       disembl_COILS  \\\n",
       "0  imtapssfeqfkvamnylqLYNVPDCLEDIQDADCSSSKCSSSASS...   \n",
       "1  TYYAWKHELLGSGTCPALPPREVLGMEELEKLPEEQvaeeelecsa...   \n",
       "2  SLHARPPQFTraqwfaiqHISLNPPRCTiamrainnyrwrcknqnt...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                      disembl_REM465  \\\n",
       "0  imtapssfeqfkvamnylqlynvpdclediqDADCSSSKCSSSASS...   \n",
       "1  tyyawkhellgsgtcpalpprevlgmeeleklpEEQVAEEELecsa...   \n",
       "2  slharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMitrpSVSQLGEGKCPSSQHL...   \n",
       "\n",
       "                                    disembl_HOTLOOPS  disembl_COILS_frac  \\\n",
       "0  IMTAPSSFeqfkvamnylqlynvpdclediqdadcSSSKCSSSASS...            0.769784   \n",
       "1  TYYAWKHELLGSGTCPALPPrevlgmeeleklpeeqvaeeelecsa...            0.791667   \n",
       "2  SLHARPPQftraqwfaiqhislnpprctiamrainnyrwrcknqnt...            0.720588   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...            0.617886   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQlgegkcpssqhl...            0.693548   \n",
       "\n",
       "         ...         disembl_HOTLOOPS_frac  ronn_avg  \\\n",
       "0        ...                      0.410072  0.495827   \n",
       "1        ...                      0.388889  0.488125   \n",
       "2        ...                      0.338235  0.368750   \n",
       "3        ...                      0.552846  0.571382   \n",
       "4        ...                      0.274194  0.556613   \n",
       "\n",
       "                                        ronn_results ronn_frac  \\\n",
       "0  imtapssfeqfkvamnylqlynvpdclediqDADCSSSKCSSSASS...  0.482014   \n",
       "1  tyyawkhellGSGTcPALPPREVLGMEELEKLPEEQVAEEELECSa...  0.493056   \n",
       "2  SLharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...  0.058824   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...  0.536585   \n",
       "4  HRKEPGARLEATRGAARPHKQGTKPMITRPSVSQLGEGKCPSSQHL...  0.620968   \n",
       "\n",
       "                                    disopred_results disopred_avg  \\\n",
       "0  ImtapssfeqfkvamnylqlynvpdclediqdadcsSSKCSSSASS...     0.387050   \n",
       "1  tyyawkhellgsgtCPALPPREVLGMEELEKLPEEQVAEEELECSA...     0.421597   \n",
       "2  slharppqftraqwfaiqhislnpprctiamrainnyrwrcknqnt...     0.034338   \n",
       "3  ARALNESKRVNNGNTAPEDSSPAKKTRRCQRQESKKMPVAGGKANK...     0.425366   \n",
       "4  HRKEPGARLEATRGAARPHKQGTkpmitrpsvsqlgegkcpssqhl...     0.236613   \n",
       "\n",
       "                                 disopred_pb_results disopred_pb_avg  \\\n",
       "0  -...................................^^^^^^^^^^...        0.317770   \n",
       "1  ..............^^^^^^--^^^^^^^^^^^^^^^^^^^-^^^^...        0.352917   \n",
       "2  .................................................        0.000000   \n",
       "3  ^^^^^^^^^^^^^^^-^---^-^^^^^^^^^^^^^--^^^^^^^^^...        0.295122   \n",
       "4  ^^^^^^^^^^^^^^^^^^^^^^^^.........................        0.185726   \n",
       "\n",
       "   disopred_frac  disopred_pb_frac  \n",
       "0       0.388489          0.359712  \n",
       "1       0.465278          0.423611  \n",
       "2       0.000000          0.000000  \n",
       "3       0.430894          0.373984  \n",
       "4       0.201613          0.209677  \n",
       "\n",
       "[5 rows x 21 columns]"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DF_disorder_features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DISOPRED3\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run DISOPRED3 on all prESTs\n",
    "\n",
    "**ANAND: (may need to fix script)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DISOPRED_DIR = MISC_DIR+'disopred_out/'\n",
    "\n",
    "disopred_cmd='/home/anand/Documents/HT_Expression/prediction_software/DISOPRED/run_disopred.pl'\n",
    "\n",
    "#subprocess.call(['bash',MISC_DIR+'disopred_out/disopred_script.sh',FASTA_DIR+'aa_sequences',disopred_cmd])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse DISOPRED disorder predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "for i,row in tqdm(DF_disorder_features.iterrows()):\n",
    "    with open(DISOPRED_DIR+str(row.prest_id)+'.diso', 'rb') as f:\n",
    "        lines = f.readlines()\n",
    "        frac = []\n",
    "        disorder_string = ''\n",
    "        for line in lines:\n",
    "            if not line.startswith('#'):\n",
    "                if line[8] == '*':\n",
    "                    disorder_string += line[6].upper()\n",
    "                else:\n",
    "                    disorder_string += line[6].lower()\n",
    "                frac.append(float(line[10:14]))\n",
    "        DF_disorder_features.loc[i,'disopred_results'] = disorder_string\n",
    "        DF_disorder_features.loc[i,'disopred_avg'] = np.mean(frac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse DISOPRED protein binding predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "for i,row in tqdm(DF_disorder_features.iterrows()):\n",
    "    with open(DISOPRED_DIR+str(row.prest_id)+'.pbdat', 'rb') as f:\n",
    "        lines = f.readlines()\n",
    "        frac = []\n",
    "        disorder_string = ''\n",
    "        for line in lines:\n",
    "            if not line.startswith('#'):\n",
    "                disorder_string += line[8]\n",
    "                try:\n",
    "                    frac.append(float(line[10:]))\n",
    "                except:\n",
    "                    frac.append(0)\n",
    "        DF_disorder_features.loc[i,'disopred_pb_results'] = disorder_string\n",
    "        DF_disorder_features.loc[i,'disopred_pb_avg'] = np.mean(frac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate numerical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_disorder_features.loc[:,'disopred_frac'] = [np.true_divide(count_upper(seq),len(seq)) for seq in DF_disorder_features.disopred_results]\n",
    "DF_disorder_features.loc[:,'disopred_pb_frac'] = [np.true_divide(seq.count('^'),len(seq)) for seq in DF_disorder_features.disopred_pb_results]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_disorder_features.to_csv(DATA_DIR+'protein_disorder.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compile Features\n",
    "[Back to Top](#Table-of-Contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DF_prest = pd.read_csv(DATA_DIR+'DF_prest.csv',index_col=0)\n",
    "DF_mRNA_features = pd.read_csv(DATA_DIR+'mRNA_features.csv',index_col=0)\n",
    "DF_physical_features = pd.read_csv(DATA_DIR+'physical_features.csv',index_col=0)\n",
    "DF_ss_features = pd.read_csv(DATA_DIR+'secondary_structure.csv',index_col=0)\n",
    "DF_disorder_features = pd.read_csv(DATA_DIR+'protein_disorder.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DF_prest_features = DF_prest.merge(DF_mRNA_features).merge(DF_physical_features).merge(DF_ss_features).merge(DF_disorder_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "split_size = len(DF_prest_features)/4\n",
    "\n",
    "DF_prest_features[:split_size].to_csv(DATA_DIR+'DF_prest_features_1.csv')\n",
    "DF_prest_features[split_size:split_size*2].to_csv(DATA_DIR+'DF_prest_features_2.csv')\n",
    "DF_prest_features[split_size*2:split_size*3].to_csv(DATA_DIR+'DF_prest_features_3.csv')\n",
    "DF_prest_features[split_size*3:].to_csv(DATA_DIR+'DF_prest_features_4.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}