{ "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 }