nuc_dynamics_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example Annealing Protocol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below function illustrates the basic operation of particle simulated annealing by implementing a simple protocol. This operates on input distrance restraints, which are specified in terms of the indices of the restrained particles and the corresponding distance limits for those restraints. The function uses the <tt>runDynamics()</tt> function to update the particle coordinates. - This function comes from the <tt>dyn_util</tt> module (written in Cython), which is assumed to have been compiled (see the <tt>setup_cython.py</tt> script) and be accessible in the same directory as this script, or otherwise on the PYTHONPATH. For full genome structure calculation from Hi-C derived restraints the <tt>nuc_dynamics.anneal_genome()</tt> should be used as this both makes the distance restraints and implements an annealing protocol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import dyn_util\n",
    "\n",
    "def anneal_structure(coords, restraintIndices, restraintLimits,\n",
    "                    tempMax=5000.0, tempMin=10.0, tempSteps=500,\n",
    "                    dynSteps=100, timeStep=0.001, randSeed=None):\n",
    "  \"\"\"\n",
    "  A simple single-stage annealing protocol with exponential\n",
    "  temperature decay and sigmoidal repulsion scheme.\n",
    "  \"\"\"\n",
    "    \n",
    "  from math import atan, log, exp, pi\n",
    "    \n",
    "  # Consider masses and fixed particles\n",
    "      \n",
    "  if randSeed is None:\n",
    "    np.random.seed()\n",
    "  else:\n",
    "    np.random.seed(randSeed)\n",
    "  \n",
    "  # Ensure inputs are numpy.ndarray  \n",
    "  coords = np.array(coords)\n",
    "  restraintIndices = np.array(restraintIndices, dtype=np.int32)\n",
    "  restraintLimits = np.array(restraintLimits)\n",
    "    \n",
    "  nPoints = len(coords)\n",
    "  \n",
    "  # Masses and radii not used, just set to 1.0\n",
    "  masses = radii = np.ones(nPoints, float)\n",
    "    \n",
    "  # Ambiguiuty not used, but represents the stride of sequential, grouped restraints\n",
    "  ambiguity = np.ones(len(restraintIndices), np.int32)\n",
    "\n",
    "  adj = 1.0 / atan(10.0)\n",
    "  decay = log(tempMax/tempMin)    \n",
    "\n",
    "  for step in range(tempSteps):\n",
    "    frac = step/float(tempSteps)\n",
    "    \n",
    "    # exponential temp decay\n",
    "    temp = tempMax * exp(-decay*frac)\n",
    "    \n",
    "    # sigmoidal repusion scheme\n",
    "    repulse = 0.5 + adj * atan(frac*20.0-10) / pi \n",
    "    \n",
    "    # update coordinates\n",
    "    dyn_util.runDynamics(coords, masses, radii, restraintIndices, restraintLimits,\n",
    "                         ambiguity, temp, timeStep, dynSteps, repulse,\n",
    "                         repDist=1.5, printInterval=np.int32(100))\n",
    "  \n",
    "    # Center\n",
    "    coords -= coords.mean(axis=0)\n",
    "  \n",
    "  return coords\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below code tests the above function by performing a small, quick calculation and demonstrates how restraints are passed into the\n",
    "function as particle indices and restraint distances.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import random\n",
    "\n",
    "# Random particle coordinates - 20 test points of (x,y,z) \n",
    "nCoords = 20\n",
    "coords = random.uniform(-1.0, 1.0, (nCoords,3)) \n",
    "  \n",
    "# Each pair represents a connection between a pair of particles\n",
    "# Restraint indices refer to indices in coordinate array\n",
    "rIndices = [(0, 1), (0, 5), (0, 18), (1, 0), (1, 2),\n",
    "            (1, 9), (2, 1), (2, 3), (2, 10), (3, 2),\n",
    "            (3, 4), (3, 17), (4, 3), (4, 5), (4, 11),\n",
    "            (5, 0), (5, 4), (5, 12), (6, 7), (6, 17),\n",
    "            (6, 19), (7, 6), (7, 14), (7, 15), (7, 16),\n",
    "            (8, 18), (9, 1), (10, 2), (11, 4), (12, 5),\n",
    "            (13, 17), (14, 7), (15, 7), (16, 7), (17, 3),\n",
    "            (17, 6), (17, 13), (18, 0), (18, 8), (19, 6)]\n",
    "  \n",
    "# Distances - upper and lower bounds for each particle\n",
    "# Unform upper and lower bounds for test restraints\n",
    "lower = 0.8\n",
    "upper = 1.2\n",
    "rDists = [(0.8, 1.2)] * len(rIndices)\n",
    "\n",
    "# Run calculation\n",
    "coords = anneal_structure(coords, rIndices, rDists)\n",
    "  \n",
    "# Show results\n",
    "for i, c in enumerate(coords):\n",
    "  print('%3d' % i, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genome Structure Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although whole genome structure calculations for Hi-C contact data may be convieniently calculated using the <tt>nuc_dynamics</tt> command-line tool, all of the functionality is also available in the <tt>nuc_dynamics.py</tt> Python module. The below example illustrates how to access the functions from this module to run a genome calculation from a Python script, though we don't recommend using this Jupyter notebook to run long calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from time import time\n",
    "\n",
    "def demo_calc_genome_structure():\n",
    "  \"\"\"\n",
    "  Example of settings for a typical genome structure calculation from input single-cell\n",
    "  Hi-C contacts stored in an NCC format file (as output from NucProcess)\n",
    "  \"\"\"\n",
    "  \n",
    "  # Import the nuc_dynamics module, which whold be in the same directory as this script\n",
    "  # or otherwise on the PYTHONPATH\n",
    "  from nuc_dynamics import calc_genome_structure\n",
    "  \n",
    "  # Set the input NCC file path (containing contact data) and the output PDB format file \n",
    "  ncc_file_path = 'example_chromo_data/Cell_1_contacts.ncc'\n",
    "  pdb_file_path = 'example_chromo_data/Cell_1_structure.pdb'\n",
    "  \n",
    "  # Number of alternative conformations to generate from repeat calculations\n",
    "  # with different random starting coordinates\n",
    "  num_models = 2\n",
    "\n",
    "  # Parameters to setup restraints and starting coords\n",
    "  general_calc_params = {'dist_power_law':-0.33,\n",
    "                         'contact_dist_lower':0.8, 'contact_dist_upper':1.2,\n",
    "                         'backbone_dist_lower':0.1, 'backbone_dist_upper':1.1,\n",
    "                         'random_seed':int(time()), 'random_radius':10.0}\n",
    "\n",
    "  # Annealing & dyamics parameters: the same for all stages\n",
    "  # (this is cautious, but not an absolute requirement)\n",
    "  anneal_params = {'temp_start':5000.0, 'temp_end':10.0, 'temp_steps':500,\n",
    "                   'dynamics_steps':100, 'time_step':0.001}\n",
    "\n",
    "  # Hierarchical scale protocol: calculations will initially use 8 Mb particles\n",
    "  # but deminish to 100 kb at the end. The whole annealing protocol (hot to cold)\n",
    "  # will be run at each size, but subsequent stages will start from the previous\n",
    "  # structure\n",
    "  particle_sizes = [8e6, 4e6, 2e6, 4e5, 2e5, 1e5]\n",
    "  \n",
    "  # Contacts must be clustered with another within this separation threshold\n",
    "  # (at both ends) to be considered supported, i.e. not isolated\n",
    "  # This removes noise contacts\n",
    "  isolation_threshold=2e6\n",
    "  \n",
    "  # Actually run the calculation with the specified parameters and input\n",
    "  # The below function will automatically create the appropriate distance restraints\n",
    "  calc_genome_structure(ncc_file_path, pdb_file_path, general_calc_params, anneal_params,\n",
    "                        particle_sizes, num_models, isolation_threshold)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test function can be run with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "demo_calc_genome_structure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}