HMM_DNA_motifs_Viterbi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ATCGG\n",
    "#formulating problem - Two states (DISEASE and NORMAL) - Controlled by SNPs at a particular motif(transcript) say ATCGG. \n",
    "#Considering background with no bias.\n",
    "#following must be calculated -\n",
    "#Probability of individuals in disease state and normal state (inititation probabilities)\n",
    "#Emission probabilities of A,T,C,G in either state.\n",
    "#transition probabilities from D-D,D-N,N-N,N-D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         0          1          2          3\n",
      "disease: 0.180000 0.018000 0.003600 0.000180\n",
      "normal: 0.100000 0.012500 0.001563 0.000195\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['normal', 'disease', 'disease', 'disease']"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#First we define all states, bases and probabilities (initiation, transition and emission)\n",
    "states = ['disease', 'normal']\n",
    "trans_prob = {'disease': {'disease':0.5, 'normal':0.5}, 'normal':{'normal':0.5,'disease':0.5}}\n",
    "init_prob = {'disease': 0.6, 'normal':0.4}\n",
    "emit_prob = {'disease':{'A':0.3,'T':0.2,'C':0.4,'G':0.1},\n",
    "    'normal':{'A':0.25,'T':0.25,'C':0.25,'G':0.25}}\n",
    "base = ('A', 'T', 'C','G')\n",
    "#creating a matrix to to store transition probabilities\n",
    "Matrix=[{}]\n",
    "#function to calculate maximum transition probabilities\n",
    "def funcalgo(states, trans_prob,init_prob,emit_prob,base):\n",
    "    for s in states:\n",
    "        Matrix[0][s]={'prob':init_prob[s]*emit_prob[s][base[0]], 'prev':None}\n",
    "        #initiation probability\n",
    "    for numobs in range(1,len(base)):\n",
    "        Matrix.append({})\n",
    "        for st in states:\n",
    "            max_prob_trans=Matrix[numobs-1][st]['prob']*trans_prob[states[0]][st]\n",
    "            prev_state_opt=states[0]\n",
    "            #keep next state fixed and vary previous state and find the max prob\n",
    "            for prev_state in states[1:]:\n",
    "                       max_prob=Matrix[numobs-1][prev_state]['prob']*trans_prob[prev_state][st]\n",
    "                       if max_prob > max_prob_trans:\n",
    "                           max_prob_trans=max_prob\n",
    "                           prev_state_opt=prev_state\n",
    "            \n",
    "            max_prob_trans_emit=max_prob_trans*emit_prob[st][base[numobs]]\n",
    "            Matrix[numobs][st]={'prob':max_prob_trans_emit,'prev':prev_state_opt}\n",
    "            \n",
    "    for line in probtable(Matrix):\n",
    "            print(line)\n",
    "\n",
    "#printing probabilities with states\n",
    "def probtable(Matrix):\n",
    "            print( \" \".join((\"%10d\" % i) for i in range(len(Matrix))))\n",
    "            for state in Matrix[0]:\n",
    "                yield \"%s: \" % state + \" \".join(\"%s\" % (\"%f\" % v[state][\"prob\"]) for v in Matrix)\n",
    "                        \n",
    "#call to algo                        \n",
    "funcalgo(states, trans_prob,init_prob,emit_prob,base)\n",
    "\n",
    "#find the state with max prob in last observable base\n",
    "last_max_prob=max(d['prob'] for d in Matrix[-1].values())\n",
    "opt=[]\n",
    "for st,value in Matrix[-1].items():\n",
    "    if value['prob']==last_max_prob:\n",
    "        opt.append(st)\n",
    "        previous=st\n",
    "        break\n",
    "\n",
    "Matrix[-1].items()\n",
    "\n",
    "#backtracking with previous states\n",
    "for t in range(len(base)-2,-1,-1):\n",
    "    opt.append(Matrix[t+1][previous]['prev'])\n",
    "    previous=Matrix[t+1][previous]['prev']\n",
    "\n",
    "#printing optimal parse\n",
    "opt\n",
    "#for t in range(0,len(base)+1):\n",
    "    #if t==0:\n",
    "        #print(\"the optimal parse of sequence is\")\n",
    "    #else:\n",
    "      #  print(opt[t-1]+\" \"+''.join(\"-\"),end='')\n",
    "    \n",
    "    \n",
    "\n",
    "        \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['normal', 'disease', 'disease', 'disease']\n",
      "           0            1            2            3\n",
      "disease: 0.18000 0.01800 0.00450 0.00022\n",
      "normal: 0.10000 0.02250 0.00281 0.00056\n",
      "dict_items([('disease', {'prob': 0.000225, 'prev': 'disease'}), ('normal', {'prob': 0.0005625, 'prev': 'disease'})])\n"
     ]
    }
   ],
   "source": [
    "obs = ('A', 'T', 'C','G')\n",
    "states = ('disease', 'normal')\n",
    "start_p = {'disease': 0.6, 'normal':0.4}\n",
    "trans_p = {'disease': {'disease':0.5, 'normal':0.5}, 'normal':{'normal':0.5,'disease':0.5}\n",
    "   }\n",
    "emit_p = {'disease':{'A':0.3,'T':0.2,'C':0.4,'G':0.1},\n",
    "    'normal':{'A':0.25,'T':0.25,'C':0.25,'G':0.25}\n",
    "   }\n",
    "\n",
    "def viterbi(obs, states, start_p, trans_p, emit_p):\n",
    "    V = [{}]\n",
    "    for st in states:\n",
    "        V[0][st] = {\"prob\": start_p[st] * emit_p[st][obs[0]], \"prev\": None}\n",
    "    # Run Viterbi when t > 0\n",
    "    for t in range(1, len(obs)):\n",
    "        V.append({})\n",
    "        for st in states:\n",
    "            max_tr_prob = V[t-1][states[0]][\"prob\"]*trans_p[states[0]][st]\n",
    "            prev_st_selected = states[0]\n",
    "            for prev_st in states[1:]:\n",
    "                tr_prob = V[t-1][prev_st][\"prob\"]*trans_p[prev_st][st]\n",
    "                if tr_prob > max_tr_prob:\n",
    "                    max_tr_prob = tr_prob\n",
    "                    prev_st_selected = prev_st\n",
    "                    \n",
    "            max_prob = max_tr_prob * emit_p[st][obs[t]]\n",
    "            V[t][st] = {\"prob\": max_prob, \"prev\": prev_st_selected}\n",
    "                    \n",
    "    for line in dptable(V):\n",
    "        print (line)\n",
    "    opt = []\n",
    "    # The highest probability\n",
    "    max_prob = max(value[\"prob\"] for value in V[-1].values())\n",
    "    previous = None\n",
    "    # Get most probable state and its backtrack \n",
    "    print(V[3].items())\n",
    "    for st, data in V[-1].items():\n",
    "        if data[\"prob\"] == max_prob:\n",
    "            opt.append(st)\n",
    "            previous = st\n",
    "            break\n",
    "    # Follow the backtrack till the first observation\n",
    "    for t in range(len(V) - 2, -1, -1):\n",
    "        opt.insert(0, V[t + 1][previous][\"prev\"])\n",
    "        previous = V[t + 1][previous][\"prev\"]\n",
    "        \n",
    "\n",
    "print(opt)\n",
    "\n",
    "def dptable(V):\n",
    "    # Print a table of steps from dictionary\n",
    "    yield \" \".join((\"%12d\" % i) for i in range(len(V)))\n",
    "    for state in V[0]:\n",
    "        yield \"%.7s: \" % state + \" \".join(\"%.7s\" % (\"%f\" % v[state][\"prob\"]) for v in V)\n",
    "        \n",
    "viterbi(obs,\n",
    "        states,\n",
    "        start_p,\n",
    "        trans_p,\n",
    "        emit_p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}