{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import zeros\n",
    "\n",
    "def exampleCost(xc, yc):\n",
    "    \"\"\" Cost function assigning 0 to match, 2 to transition, 4 to\n",
    "        transversion, and 8 to a gap \"\"\"\n",
    "    if xc == yc: return 0 # match\n",
    "    if xc == '-' or yc == '-': return 8 # gap\n",
    "    minc, maxc = min(xc, yc), max(xc, yc)\n",
    "    if minc == 'A' and maxc == 'G': return 2 # transition\n",
    "    elif minc == 'C' and maxc == 'T': return 2 # transition\n",
    "    return 4 # transversion\n",
    "\n",
    "def globalAlignment(x, y, s):\n",
    "    \"\"\" Calculate global alignment value of sequences x and y using\n",
    "        dynamic programming.  Return global alignment value. \"\"\"\n",
    "    D = zeros((len(x)+1, len(y)+1), dtype=int)\n",
    "    for j in range(1, len(y)+1):\n",
    "        D[0, j] = D[0, j-1] + s('-', y[j-1])\n",
    "    for i in range(1, len(x)+1):\n",
    "        D[i, 0] = D[i-1, 0] + s(x[i-1], '-')\n",
    "    for i in range(1, len(x)+1):\n",
    "        for j in range(1, len(y)+1):\n",
    "            D[i, j] = min(D[i-1, j-1] + s(x[i-1], y[j-1]), # diagonal\n",
    "                          D[i-1, j  ] + s(x[i-1], '-'),    # vertical\n",
    "                          D[i  , j-1] + s('-',    y[j-1])) # horizontal\n",
    "    return D, D[len(x), len(y)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D, globalAlignmentValue = globalAlignment('TACGTCAGC', 'TATGTCATGC', exampleCost)\n",
    "globalAlignmentValue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  8 16 24 32 40 48 56 64 72 80]\n",
      " [ 8  0  8 16 24 32 40 48 56 64 72]\n",
      " [16  8  0  8 16 24 32 40 48 56 64]\n",
      " [24 16  8  2 10 18 24 32 40 48 56]\n",
      " [32 24 16 10  2 10 18 26 34 40 48]\n",
      " [40 32 24 16 10  2 10 18 26 34 42]\n",
      " [48 40 32 24 18 10  2 10 18 26 34]\n",
      " [56 48 40 32 26 18 10  2 10 18 26]\n",
      " [64 56 48 40 32 26 18 10  6 10 18]\n",
      " [72 64 56 48 40 34 26 18 12 10 10]]\n"
     ]
    }
   ],
   "source": [
    "print(D)"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}