variant_repeat_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10, 10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deletion_score = -1\n",
    "gap_score = -1\n",
    "mismatch_score = -1\n",
    "match_score = +1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def align(template, sequence, cache = None):\n",
    "    if cache is None:\n",
    "        cache = dict()\n",
    "    if (template, sequence) in cache:\n",
    "        return cache[(template, sequence)]\n",
    "    if not template:\n",
    "        return len(sequence) * deletion_score, \"\"\n",
    "    if not sequence:\n",
    "        return len(template) * gap_score, template\n",
    "    \n",
    "    template_head, template_tail = template[0], template[1:]\n",
    "    sequence_head, sequence_tail = sequence[0], sequence[1:]\n",
    "    \n",
    "    if template_head == sequence_head:\n",
    "        match_mismatch_score = match_score\n",
    "    else:\n",
    "        match_mismatch_score = mismatch_score\n",
    "    \n",
    "    s, seq = align(template_tail, sequence_tail, cache)\n",
    "    match_result = s + match_mismatch_score, template_head + seq\n",
    "    s, seq = align(template, sequence_tail, cache)\n",
    "    delete_result = s + deletion_score, seq\n",
    "    s, seq = align(template_tail, sequence, cache)\n",
    "    gap_result = s + gap_score, template_head + seq\n",
    "    \n",
    "    results = [match_result, delete_result, gap_result]\n",
    "    results.sort()\n",
    "    result = results[-1]\n",
    "    cache[(template, sequence)] = result\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 'GCATGCT')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "align(\"GCATGCT\", \"GATTACA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def produce_sequences(pattern, max_length, outer_index, inner_index, partial_result):\n",
    "    if len(partial_result) >= max_length:\n",
    "        return [partial_result]\n",
    "\n",
    "    current_pattern = pattern[outer_index]\n",
    "    results = []\n",
    "    next_pattern_available = outer_index != len(pattern) - 1 and (inner_index == 0)\n",
    "        \n",
    "    next_inner_index = (inner_index + 1) % len(current_pattern)\n",
    "    current_character = pattern[outer_index][inner_index]\n",
    "    next_partial = partial_result + current_character\n",
    "    \n",
    "    if next_pattern_available:\n",
    "        jump_right_now = produce_sequences(pattern, max_length, outer_index + 1, 0, partial_result)\n",
    "        for p in jump_right_now:\n",
    "            results.append(p)\n",
    "    \n",
    "    same_pattern_result = produce_sequences(pattern, max_length, outer_index, next_inner_index, next_partial)\n",
    "    for p in same_pattern_result:\n",
    "        results.append(p)\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def produce_sequences_outer(pattern, max_length):\n",
    "    return produce_sequences(pattern, max_length, 0, 0, \"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BBBBB', 'AABBB', 'AAAAB', 'AAAAA']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "produce_sequences_outer((\"AA\", \"BB\"), 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['CTGCTGCTGCTG',\n",
       " 'CCGCTGCTGCTG',\n",
       " 'CCGCTGCCGCTG',\n",
       " 'CTGCTGCTGCTG',\n",
       " 'CTGCCGCTGCTG',\n",
       " 'CTGCCGCTGCCG',\n",
       " 'CTGCTGCTGCTG',\n",
       " 'CTGCTGCCGCTG',\n",
       " 'CTGCTGCTGCTG',\n",
       " 'CTGCTGCTGCCG',\n",
       " 'CTGCTGCTGCTG']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "produce_sequences_outer((\"CTG\", \"CCGCTG\", \"CTG\"), 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def produce_sequence_and_align(pattern, template):\n",
    "    sequences = produce_sequences_outer(pattern, len(template))\n",
    "    results = sorted([align(seq, template) for seq in sequences])\n",
    "    return results[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def collapse_sequence(pattern, sequence):\n",
    "    answer = []\n",
    "    repetitions = 0\n",
    "    for elem in pattern:\n",
    "        while True:\n",
    "            tried = elem * repetitions\n",
    "            if sequence[:len(tried)] != tried:\n",
    "                answer.append(repetitions - 1)\n",
    "                sequence = sequence[len(tried) - len(elem):]\n",
    "                repetitions = 0\n",
    "                break\n",
    "            repetitions += 1\n",
    "    return list(zip(answer, pattern))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def best_templates_from_raw_reads(pattern, sequences):\n",
    "    avg_sequence_length = round(sum(len(sequence) for sequence in sequences)/len(sequences))\n",
    "    candidate_patterns = produce_sequences_outer(pattern, avg_sequence_length)\n",
    "    scores = {sequence: 0 for sequence in candidate_patterns}\n",
    "    for pattern in candidate_patterns:\n",
    "        for sequence in sequences:\n",
    "            scores[pattern] += align(pattern, sequence)[0]\n",
    "    \n",
    "    for pattern in candidate_patterns:\n",
    "        for sequence in sequences:\n",
    "            scores[pattern] += align(pattern, sequence)[0]\n",
    "    \n",
    "    sorted_sequences = sorted([(value, sequence) for sequence, value in scores.items()], reverse=True)\n",
    "    return sorted_sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo on non-biological data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pattern = [\"AA\", \"ABA\", \"BB\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BBBBBBB',\n",
       " 'ABABBBB',\n",
       " 'ABAABAB',\n",
       " 'ABAABAA',\n",
       " 'AABBBBB',\n",
       " 'AAABABB',\n",
       " 'AAABAAB',\n",
       " 'AAAABBB',\n",
       " 'AAAAABA',\n",
       " 'AAAAAAB',\n",
       " 'AAAAAAA',\n",
       " 'AAAAAAA']"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "produce_sequences_outer(pattern, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(54, 'AAABABB'),\n",
       " (38, 'AAABAAB'),\n",
       " (34, 'AABBBBB'),\n",
       " (32, 'AAAABBB'),\n",
       " (30, 'ABABBBB'),\n",
       " (30, 'ABAABAB'),\n",
       " (18, 'AAAAABA'),\n",
       " (18, 'AAAAAAB'),\n",
       " (10, 'ABAABAA'),\n",
       " (-4, 'AAAAAAA'),\n",
       " (-6, 'BBBBBBB')]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reads = [\"AAABABB\", \"AACBABB\", \"AABBABB\", \"AAABBBB\", \"AAABAB\"]\n",
    "result = best_templates_from_raw_reads(pattern, reads)\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "best_template = result[0][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'AAABABB'"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 'AA'), (1, 'ABA'), (1, 'BB')]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collapse_sequence(pattern, best_template)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problems with biological data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hypothetical_pattern = [\"CTG\", \"CCGCTG\", \"CTG\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sequence = 'GCTGTTGCTGCTGCTTGCTGCTGCTTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCATGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCGCTGCTGCTGCTGCTGCTGCTGCTGCCGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGATGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCCGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCGCTGCCGCTGCCGCTGCGCTGCCGCTGCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCGCTGCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCGCTGCGCTGCGCTGCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCGCTGCCGCTGCGCTGCGCTGCGCTGCCGCTGCGCTGCCGCTGCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCGCTGCGCTGCCGCTGCCGCTGCGCTGCGCTGCCGCTGCGCTGCCGCTGCCGCTGCGCTGCCGCTGCCGCAGCCGCTGCCGCTGCCGCTGCCGCTGCCGCTGCCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCCGCTGCTGCTGCT'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1389"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "start = time.time()\n",
    "candidate_templates = produce_sequences_outer([\"CTG\", \"CTGCCG\", \"CTG\"], len(sequence))\n",
    "stop = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "took 29.27958345413208s\n"
     ]
    }
   ],
   "source": [
    "print(\"took {}s\".format(stop - start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54288"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(candidate_templates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "align(candidate_templates[0], sequence)\n",
    "stop = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "single_alignment = stop - start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "took 19.319965600967407s\n"
     ]
    }
   ],
   "source": [
    "print(\"took {}s\".format(single_alignment))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_time_seconds = single_alignment * len(candidate_templates) * 700"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "734189604.781723"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_time_seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_years = total_time_seconds/(3600*24*30*12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23.604346861552308"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_years"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Biostatisticians in Glasgow suggested I use a \"Markov Change-Point model\" and build a maximum liklihood model (should be faster).\n",
    "2. I haven't tried it yet, but interested to collaborate!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 Clinical Applicatoins",
   "language": "python",
   "name": "clinical_applications"
  },
  "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": 2
}