Quick check of contact distances for NtrC

ntrc_test.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:61cb79ce5f2705b13a2aad4c3913d9f6ef09be5eae0cbe5a1359f6d05ac55941"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import mdtraj as md"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# fetch pdbs from RCSB\n",
      "!wget http://www.rcsb.org/pdb/files/1DC7.pdb -O 1DC7.pdb\n",
      "!wget http://www.rcsb.org/pdb/files/1DC8.pdb -O 1DC8.pdb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "--2015-04-07 11:02:00--  http://www.rcsb.org/pdb/files/1DC7.pdb\r\n",
        "Resolving www.rcsb.org... 128.6.70.10\r\n",
        "Connecting to www.rcsb.org|128.6.70.10|:80... connected.\r\n",
        "HTTP request sent, awaiting response... "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "200 OK\r\n",
        "Length: unspecified [text/plain]\r\n",
        "Saving to: \u20181DC7.pdb\u2019\r\n",
        "\r\n",
        "\r",
        "    [<=>                                    ] 0           --.-K/s              "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "    [ <=>                                   ] 168,642     --.-K/s   in 0.05s   \r\n",
        "\r\n",
        "2015-04-07 11:02:00 (2.93 MB/s) - \u20181DC7.pdb\u2019 saved [168642]\r\n",
        "\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "--2015-04-07 11:02:00--  http://www.rcsb.org/pdb/files/1DC8.pdb\r\n",
        "Resolving www.rcsb.org... 128.6.70.10\r\n",
        "Connecting to www.rcsb.org|128.6.70.10|:80... connected.\r\n",
        "HTTP request sent, awaiting response... "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "200 OK\r\n",
        "Length: unspecified [text/plain]\r\n",
        "Saving to: \u20181DC8.pdb\u2019\r\n",
        "\r\n",
        "\r",
        "    [<=>                                    ] 0           --.-K/s              "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "    [ <=>                                   ] 172,449     --.-K/s   in 0.06s   \r\n",
        "\r\n",
        "2015-04-07 11:02:00 (2.62 MB/s) - \u20181DC8.pdb\u2019 saved [172449]\r\n",
        "\r\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def load_pdb(filename):\n",
      "    '''Load pdb file and extract CA atoms'''\n",
      "    traj = md.load(filename)\n",
      "    ca_indices = traj.top.select('name CA')\n",
      "    \n",
      "    return traj.atom_slice(ca_indices)\n",
      "\n",
      "def calc_distances(traj):\n",
      "    '''Calc pairwise distances and return in angstroms in both square form and unique atom pairs\n",
      "    as well as the atom index pairs'''\n",
      "    atom_pairs = []\n",
      "    for i in xrange(traj.n_atoms - 1):\n",
      "        for j in xrange(i+1, traj.n_atoms):\n",
      "            atom_pairs.append((i,j))\n",
      "    \n",
      "    d = md.compute_distances(traj, atom_pairs=atom_pairs, periodic=False)\n",
      "    \n",
      "    # Convert to square matrix\n",
      "    contact_dist = md.geometry.squareform(d, atom_pairs)\n",
      "    \n",
      "    return 10*contact_dist[0], 10*d[0], atom_pairs\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load pdb file\n",
      "stateA = load_pdb('1DC7.pdb')\n",
      "stateB = load_pdb('1DC8.pdb')\n",
      "\n",
      "assert stateA.n_atoms == 124\n",
      "assert stateB.n_atoms == 124"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "distA, dA, atom_pairs = calc_distances(stateA)\n",
      "distB, dB, atom_pairs = calc_distances(stateB)\n",
      "\n",
      "assert len(atom_pairs) == int(0.5*124*123)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the number of CA-CA pairs with distances less than 11.5 angstroms\n",
      "np.sum(distA < 11.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "3184"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(distB < 11.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "3242"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the number of CA-CA pairs with distances less than 11.5 angstroms in both states\n",
      "np.sum(np.logical_and((distA < 11.5), (distB < 11.5)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "2626"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the number of CA-CA pairs with distances less than 11.5 angstroms in either states\n",
      "np.sum(np.logical_or((distA < 11.5), (distB < 11.5)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "3800"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now just look at unique atom pairs since the above double counts and include $i,j$ and $j,i$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(dA < 11.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "1530"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.sum(dB < 11.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "1559"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the number of CA-CA pairs with distances less than 11.5 angstroms in both states\n",
      "np.sum(np.logical_and((dA < 11.5), (dB < 11.5)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "1251"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the number of CA-CA pairs with distances less than 11.5 angstroms in either states\n",
      "np.sum(np.logical_or((dA < 11.5), (dB < 11.5)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "1838"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Number of pairs whose delta dist is greater than 5 angstroms\n",
      "np.sum(np.abs(dA - dB) > 5.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "553"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Number of pairs whose delta dist is greater than 5 angstroms and whose d_ij is less than 11.5 in both states\n",
      "indx = np.logical_and((dA < 11.5), (dB < 11.5))\n",
      "np.sum(np.abs(dA[indx] - dB[indx]) > 5.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "10"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Number of pairs whose delta dist is greater than 5 angstroms and whose d_ij is less than 11.5 in either states\n",
      "indx = np.logical_or((dA < 11.5), (dB < 11.5))\n",
      "np.sum(np.abs(dA[indx] - dB[indx]) > 5.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "126"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}