2.03_ApproximateMatching.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import string\n",
    "\n",
    "def z_array(s):\n",
    "    \"\"\" Use Z algorithm (Gusfield theorem 1.4.1) to preprocess s \"\"\"\n",
    "    assert len(s) > 1\n",
    "    z = [len(s)] + [0] * (len(s)-1)\n",
    "\n",
    "    # Initial comparison of s[1:] with prefix\n",
    "    for i in range(1, len(s)):\n",
    "        if s[i] == s[i-1]:\n",
    "            z[1] += 1\n",
    "        else:\n",
    "            break\n",
    "    \n",
    "    r, l = 0, 0\n",
    "    if z[1] > 0:\n",
    "        r, l = z[1], 1\n",
    "    \n",
    "    for k in range(2, len(s)):\n",
    "        assert z[k] == 0\n",
    "        if k > r:\n",
    "            # Case 1\n",
    "            for i in range(k, len(s)):\n",
    "                if s[i] == s[i-k]:\n",
    "                    z[k] += 1\n",
    "                else:\n",
    "                    break\n",
    "            r, l = k + z[k] - 1, k\n",
    "        else:\n",
    "            # Case 2\n",
    "            # Calculate length of beta\n",
    "            nbeta = r - k + 1\n",
    "            zkp = z[k - l]\n",
    "            if nbeta > zkp:\n",
    "                # Case 2a: Zkp wins\n",
    "                z[k] = zkp\n",
    "            else:\n",
    "                # Case 2b: Compare characters just past r\n",
    "                nmatch = 0\n",
    "                for i in range(r+1, len(s)):\n",
    "                    if s[i] == s[i - k]:\n",
    "                        nmatch += 1\n",
    "                    else:\n",
    "                        break\n",
    "                l, r = k, r + nmatch\n",
    "                z[k] = r - k + 1\n",
    "    return z\n",
    "\n",
    "\n",
    "def n_array(s):\n",
    "    \"\"\" Compile the N array (Gusfield theorem 2.2.2) from the Z array \"\"\"\n",
    "    return z_array(s[::-1])[::-1]\n",
    "\n",
    "\n",
    "def big_l_prime_array(p, n):\n",
    "    \"\"\" Compile L' array (Gusfield theorem 2.2.2) using p and N array.\n",
    "        L'[i] = largest index j less than n such that N[j] = |P[i:]| \"\"\"\n",
    "    lp = [0] * len(p)\n",
    "    for j in range(len(p)-1):\n",
    "        i = len(p) - n[j]\n",
    "        if i < len(p):\n",
    "            lp[i] = j + 1\n",
    "\n",
    "    return lp\n",
    "\n",
    "\n",
    "def big_l_array(p, lp):\n",
    "    \"\"\" Compile L array (Gusfield theorem 2.2.2) using p and L' array.\n",
    "        L[i] = largest index j less than n such that N[j] >= |P[i:]| \"\"\"\n",
    "    l = [0] * len(p)\n",
    "    l[1] = lp[1]\n",
    "    for i in range(2, len(p)):\n",
    "        l[i] = max(l[i-1], lp[i])\n",
    "    return l\n",
    "\n",
    "\n",
    "def small_l_prime_array(n):\n",
    "    \"\"\" Compile lp' array (Gusfield theorem 2.2.4) using N array. \"\"\"\n",
    "    small_lp = [0] * len(n)\n",
    "    for i in range(len(n)):\n",
    "        if n[i] == i+1:  # prefix matching a suffix\n",
    "            small_lp[len(n)-i-1] = i+1\n",
    "    for i in range(len(n)-2, -1, -1):  # \"smear\" them out to the left\n",
    "        if small_lp[i] == 0:\n",
    "            small_lp[i] = small_lp[i+1]\n",
    "    return small_lp\n",
    "\n",
    "\n",
    "def good_suffix_table(p):\n",
    "    \"\"\" Return tables needed to apply good suffix rule. \"\"\"\n",
    "    n = n_array(p)\n",
    "    lp = big_l_prime_array(p, n)\n",
    "    return lp, big_l_array(p, lp), small_l_prime_array(n)\n",
    "\n",
    "\n",
    "def good_suffix_mismatch(i, big_l_prime, small_l_prime):\n",
    "    \"\"\" Given a mismatch at offset i, and given L/L' and l' arrays,\n",
    "        return amount to shift as determined by good suffix rule. \"\"\"\n",
    "    length = len(big_l_prime)\n",
    "    assert i < length\n",
    "    if i == length - 1:\n",
    "        return 0\n",
    "    i += 1  # i points to leftmost matching position of P\n",
    "    if big_l_prime[i] > 0:\n",
    "        return length - big_l_prime[i]\n",
    "    return length - small_l_prime[i]\n",
    "\n",
    "\n",
    "def good_suffix_match(small_l_prime):\n",
    "    \"\"\" Given a full match of P to T, return amount to shift as\n",
    "        determined by good suffix rule. \"\"\"\n",
    "    return len(small_l_prime) - small_l_prime[1]\n",
    "\n",
    "\n",
    "def dense_bad_char_tab(p, amap):\n",
    "    \"\"\" Given pattern string and list with ordered alphabet characters, create\n",
    "        and return a dense bad character table.  Table is indexed by offset\n",
    "        then by character. \"\"\"\n",
    "    tab = []\n",
    "    nxt = [0] * len(amap)\n",
    "    for i in range(0, len(p)):\n",
    "        c = p[i]\n",
    "        assert c in amap\n",
    "        tab.append(nxt[:])\n",
    "        nxt[amap[c]] = i+1\n",
    "    return tab\n",
    "\n",
    "\n",
    "class BoyerMoore(object):\n",
    "    \"\"\" Encapsulates pattern and associated Boyer-Moore preprocessing. \"\"\"\n",
    "    \n",
    "    def __init__(self, p, alphabet='ACGT'):\n",
    "        self.p = p\n",
    "        self.alphabet = alphabet\n",
    "\n",
    "        # Create map from alphabet characters to integers\n",
    "        self.amap = {}\n",
    "        for i in range(len(self.alphabet)):\n",
    "            self.amap[self.alphabet[i]] = i\n",
    "\n",
    "        # Make bad character rule table\n",
    "        self.bad_char = dense_bad_char_tab(p, self.amap)\n",
    "\n",
    "        # Create good suffix rule table\n",
    "        _, self.big_l, self.small_l_prime = good_suffix_table(p)\n",
    "\n",
    "    def bad_character_rule(self, i, c):\n",
    "        \"\"\" Return # skips given by bad character rule at offset i \"\"\"\n",
    "        assert c in self.amap\n",
    "        ci = self.amap[c]\n",
    "        assert i > (self.bad_char[i][ci]-1)\n",
    "        return i - (self.bad_char[i][ci]-1)\n",
    "\n",
    "    def good_suffix_rule(self, i):\n",
    "        \"\"\" Given a mismatch at offset i, return amount to shift\n",
    "            as determined by (weak) good suffix rule. \"\"\"\n",
    "        length = len(self.big_l)\n",
    "        assert i < length\n",
    "        if i == length - 1:\n",
    "            return 0\n",
    "        i += 1  # i points to leftmost matching position of P\n",
    "        if self.big_l[i] > 0:\n",
    "            return length - self.big_l[i]\n",
    "        return length - self.small_l_prime[i]\n",
    "\n",
    "    def match_skip(self):\n",
    "        \"\"\" Return amount to shift in case where P matches T \"\"\"\n",
    "        return len(self.small_l_prime) - self.small_l_prime[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def boyer_moore(p, p_bm, t):\n",
    "    \"\"\" Do Boyer-Moore matching \"\"\"\n",
    "    i = 0\n",
    "    occurrences = []\n",
    "    while i < len(t) - len(p) + 1:\n",
    "        shift = 1\n",
    "        mismatched = False\n",
    "        for j in range(len(p)-1, -1, -1):\n",
    "            if p[j] != t[i+j]:\n",
    "                skip_bc = p_bm.bad_character_rule(j, t[i+j])\n",
    "                skip_gs = p_bm.good_suffix_rule(j)\n",
    "                shift = max(shift, skip_bc, skip_gs)\n",
    "                mismatched = True\n",
    "                break\n",
    "        if not mismatched:\n",
    "            occurrences.append(i)\n",
    "            skip_gs = p_bm.match_skip()\n",
    "            shift = max(shift, skip_gs)\n",
    "        i += shift\n",
    "    return occurrences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def approximate_match(p, t, n):\n",
    "    segment_length = int(round(len(p) / (n+1)))\n",
    "    all_matches = set()\n",
    "    for i in range(n+1):\n",
    "        start = i*segment_length\n",
    "        end = min((i+1)*segment_length, len(p))\n",
    "        p_bm = BoyerMoore(p[start:end], alphabet='ACGT')\n",
    "        matches = boyer_moore(p[start:end], p_bm, t)\n",
    "        \n",
    "        # Extend matching segments to see if whole p matches\n",
    "        for m in matches:\n",
    "            if m < start or m-start+len(p) > len(t):\n",
    "                continue\n",
    "            \n",
    "            mismatches = 0\n",
    "            \n",
    "            for j in range(0, start):\n",
    "                if not p[j] == t[m-start+j]:\n",
    "                    mismatches += 1\n",
    "                    if mismatches > n:\n",
    "                        break\n",
    "            for j in range(end, len(p)):\n",
    "                if not p[j] == t[m-start+j]:\n",
    "                    mismatches += 1\n",
    "                    if mismatches > n:\n",
    "                        break\n",
    "            \n",
    "            if mismatches <= n:\n",
    "                all_matches.add(m - start)\n",
    "    return list(all_matches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 5]\n"
     ]
    }
   ],
   "source": [
    "p = 'AACTTG'\n",
    "t = 'CACTTAATTTG'\n",
    "print(approximate_match(p, t, 2))"
   ]
  }
 ],
 "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
}