BA10G - Alignment to profile HMM .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BA10G Multiple sequence alignment with profile HMM\n",
    "\n",
    "**STOP and Think**: Since the HMM diagram in the figure below has 25 nodes — not including the start and end states — the Viterbi graph for the string emitted in this figure has 25 rows. How many columns does this Viterbi graph have?\n",
    "\n",
    "*the length of the sequence, + any deletions are states too??*\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src='http://bioinformaticsalgorithms.com/images/HMM/aligning_text_against_profile.png' width = 500 px>\n",
    "\n",
    "<br>\n",
    "\n",
    "<br>\n",
    "\n",
    "We are now ready to align a string Text to a multiple alignment by constructing the Viterbi graph for this string (see figure below) and solving the Decoding Problem to find the most likely hidden path.\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src =http://bioinformaticsalgorithms.com/images/HMM/AEFDFDC_path_1.png width = 500 px >\n",
    "\n",
    "\n",
    "**Figure**: The Viterbi graph for HMM(Alignment, θ) and a path in this graph (shown in purple) corresponding to the hidden path for the emitted string AEFDFDC from the alignment we have been working with. Edges between columns correspond to allowed transitions in the HMM diagram and have an implied rightward orientation. Edges entering nodes corresponding to deletion states are dashed. Emitted symbols are shown beneath each column.\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "￼STOP and Think: Find paths through the Viterbi graph corresponding to the bottom four hidden paths in the figure below. What happens?\n",
    "\n",
    "*it just doesn't work ? there aren't any edges to the end for most sequences throught the viterbi graph* \n",
    "\n",
    "<br>\n",
    "<img src = http://bioinformaticsalgorithms.com/images/HMM/alignment_HMM_paths.png width = 400px>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#### The troublesome silent states\n",
    "\n",
    "If you reached this point without any questions about the figure two steps ago, then we have successfully concealed from you that solving **the Decoding Problem for HMMs with silent states** is not as simple as it may appear: the graph we proposed is not a Viterbi graph! \n",
    "\n",
    "\n",
    "To see why not, consider the path in the figure below, which emits the same string but passes through one fewer silent deletion state, thus reducing the number of columns by one. But the *Viterbi graph is not allowed to change depending on the hidden path π, since we know nothing about the hidden path in advance*! \n",
    "  \n",
    "\n",
    "Instead, the number of **columns in the Viterbi graph must equal the length of the emitted string, a condition that is violated in the two figures we have proposed**.\n",
    "\n",
    "￼STOP and Think: How can we modify the notion of the Viterbi graph for HMMs with silent states?\n",
    "\n",
    "<br>\n",
    "\n",
    "<img src = http://bioinformaticsalgorithms.com/images/HMM/AEFDFDC_path_2.png width=450px>\n",
    "\n",
    "**Figure**: Another path through another “Viterbi graph” emitting the same string AEFDFDC.\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "More generally, **the Viterbi algorithm does not tolerate silent states other than the initial and terminal states**. In other words, this algorithm assumes that node (k, i) in the Viterbi graph describes the event “the HMM emitted symbol xi when it was in state k”. However, if k is a silent state, then the role of the node (k, i) in the Viterbi graph is poorly defined, as it is unclear how to define the weight of edges entering this node.\n",
    "\n",
    "\n",
    "Fortunately, we can fix this issue in the case of profile HMMs by defining the Viterbi graph with |States| rows and |Text| columns (see figure below). Every time the HMM moves into a deletion state, rather than crossing over to the next column of the Viterbi graph (as in the two failed figures), we will move within the same column. When the HMM moves into a match or insertion state, we will move to the next column. As a result, every column of the Viterbi graph corresponds to a single emitted symbol, even though a path can pass through more than one state in a given column.\n",
    "\n",
    "\n",
    "**Exercise Break**: Show that the vertical edge connecting (i, l) to (i, k), where k is a deletion state, should be assigned weight equal to transitionl, k.\n",
    "\n",
    "**STOP and Think**: Take another look at the graph from the previous step, reproduced below. Are there any remaining issues with it?\n",
    "\n",
    "\n",
    "<br> \n",
    "<img src = http://bioinformaticsalgorithms.com/images/HMM/profile_HMM_vertical_deletion_edges.png width = 400 px>\n",
    "<br>\n",
    "\n",
    "**Figure**: The Viterbi graph with |States| rows and |Text| columns for the profile HMM we have been considering emitting a string Text of length 7 so that edges entering deletion states are drawn downward within the same column instead of between columns as in the two previously proposed figures. The purple path corresponds to the path through the HMM emitting AEFDFDC.\n",
    "\n",
    " \n",
    " There is still a minor flaw with the graph from the previous step. If the HMM moves from the initial state into Deletion(1), then the HMM will move through the first column without emitting a symbol. We will therefore transform the initial state into a column of silent states containing the initial state and all deletion states (figure below). This way, if the HMM enters Deletion(1) from the initial state, it can move downward through deletion states before transitioning to a match or insertion state in the first column.\n",
    "\n",
    "\n",
    "<br> \n",
    "<img src =  http://bioinformaticsalgorithms.com/images/HMM/profile_HMM_Viterbi_complete.png width = 400 px>\n",
    "<br>\n",
    "\n",
    "**Figure**: The final Viterbi graph for the profile HMM emitting a string of length 7. Edges within the same column have a downward orientation; edges between columns have a rightward orientation. Once again, the purple path corresponds to the path through the HMM emitting AEFDFDC.\n",
    "\n",
    "<br>\n",
    "\n",
    "<br>\n",
    "\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "<br>\n",
    "\n",
    "<br>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "def parse_HMM(file): # with pseudocount sigma\n",
    "    lines = f.readlines()\n",
    "    text = lines[0].strip()\n",
    "    # if error: check if tab or space in this input line\n",
    "#     (theta, sigma) = (float(i) for i in lines[2].strip().split('\\t'))\n",
    "    (theta, sigma) = (float(i) for i in lines[2].strip().split(' '))\n",
    "    alpha = lines[4].strip().split()\n",
    "    alpha = dict(zip(alpha, range(len(alpha))))\n",
    "    alignment = [line.strip() for line in lines[6:]]\n",
    "    print(\"Text: \", text)\n",
    "    print(\"Alignment:\")\n",
    "    for a in alignment:\n",
    "        print('\\t', a)\n",
    "    print(\"Alpha:\", alpha)\n",
    "    print(\"Theta:\", theta)\n",
    "    print(\"Sigma:\", sigma)\n",
    "    return text, alignment, alpha, theta, sigma\n",
    "\n",
    "### algorithm for building profile HMM from an alignment|theta\n",
    "def profileHMM(alignment, alpha, theta):\n",
    "\n",
    "    def get_seed_alignment(alignment, theta, alphabet):\n",
    "        \"\"\"returns seed alignment columns for HMM\"\"\"\n",
    "        k = len(alignment[0]); a = len(alphabet.keys())\n",
    "        freq = np.zeros(shape=(a+1, k)) \n",
    "        for seq in alignment:\n",
    "            for i in range(k):\n",
    "                if seq[i] == '-':\n",
    "                    freq[a][i] += 1\n",
    "                else:\n",
    "                    freq[alphabet[seq[i]]][i] +=1\n",
    "        n = len(alignment)\n",
    "        seed = [x/n < theta for x in freq[a]]\n",
    "        return seed\n",
    "\n",
    "    def state_transition(T, prev, kind, S):\n",
    "        \"\"\" find next Match, Del, Ins, or End state in [states] \"\"\"\n",
    "        x=0\n",
    "        if S[prev][0] == 'M':\n",
    "            x=1\n",
    "        for nxt in range(prev+1+x, len(T)):\n",
    "            if S[nxt][0] == kind[0]:\n",
    "                T[prev][nxt] += 1\n",
    "                return T, nxt\n",
    "            \n",
    "    def normalize_matrices(T,E):\n",
    "        \"\"\"Normalize Transmission and Emission frequencies to sum of row for state\"\"\"\n",
    "        for state in range(len(S)):\n",
    "            if sum(T[state]) > 0:\n",
    "                T[state] =  T[state]/sum(T[state])\n",
    "            if sum(E[state]) >0:\n",
    "                E[state] =  E[state]/sum(E[state])\n",
    "        return T, E\n",
    "    \n",
    "    def pseudo_normalize(T, E):\n",
    "        \"\"\"adds pseudocounts to HMM matrices\"\"\"\n",
    "        # add pseudos and normalize E\n",
    "        for i in range(1, len(S)-1):\n",
    "            if S[i][0] != 'D':\n",
    "                E[i] += 0.01\n",
    "                E[i] = E[i]/sum(E[i])\n",
    "        # add pseudos and normalize T\n",
    "        for i in range(len(S)-1):\n",
    "            start = 3*((i+1)//3)+1\n",
    "            end = start+3\n",
    "            for j in range(start, end):\n",
    "                if j < len(S):\n",
    "                    T[i][j] += 0.01\n",
    "            T[i] = T[i]/sum(T[i])\n",
    "        return T, E\n",
    "    \n",
    "    ### Walk each sequence to count transitions, emissions\n",
    "    seed = get_seed_alignment(alignment, theta, alpha)\n",
    "    n = len(alignment); k = len(alignment[0])\n",
    "    S = ['S','I0']+list(c+str(n) for n in range(1,sum(seed)+1) for c in \"MDI\")+['E']\n",
    "    E = np.zeros(shape=(len(S), len(alpha.keys())))\n",
    "    T = np.zeros(shape=(len(S), len(S)))\n",
    "\n",
    "    for seq in alignment:\n",
    "        # 'state' is the hidden state col/row; i is pos in alignment\n",
    "        state = 0; i = 0\n",
    "        while i < k:\n",
    "            # seed: either match or del as hidden state\n",
    "            if seed[i]:\n",
    "                if seq[i] in alpha:\n",
    "                    T, state = state_transition(T, state, 'Match', S)\n",
    "                    E[state][alpha[seq[i]]] +=1\n",
    "                else:\n",
    "                    T, state = state_transition(T, state, 'Deletion', S)\n",
    "            # not seed: either insert or nothing\n",
    "            else:\n",
    "                # count any emissions before next seed column\n",
    "                emits = []\n",
    "                while not seed[i]:\n",
    "                    if seq[i] in alpha:\n",
    "                        emits.append(seq[i])\n",
    "                    i += 1\n",
    "                    if i == k: \n",
    "                        break  # hit end of sequence\n",
    "                i-=1                \n",
    "                if len(emits)>0:\n",
    "                    # count the transition(state, insertion)\n",
    "                    T, state = state_transition(T, state, 'Insert', S)\n",
    "                    # get all emissions(state) in insertion\n",
    "                    for symbol in emits:\n",
    "                        E[state][alpha[symbol]] += 1\n",
    "                    # count all symbols as cyclic transition t(ins_x, ins_x)\n",
    "                    if len(emits) >1:\n",
    "                        T[state][state] += len(emits)-1\n",
    "                # else do nothing, just gaps in align\n",
    "            i += 1\n",
    "        # from last state to 'End'    \n",
    "        T, state = state_transition(T, state, 'End', S)\n",
    "\n",
    "    T,E = normalize_matrices(T,E)\n",
    "    T, E = pseudo_normalize(T, E)\n",
    "    return S, T, E, seed\n",
    "\n",
    "\n",
    "\n",
    "def Viterbi_alignment(T,E,S,alpha,text):\n",
    "    \"\"\"Aligns text to HMM using Viterbi algorithm for optimal hidden path\"\"\"\n",
    "    \n",
    "    # initialize the Viterbi graph (|states|*|text|)\n",
    "    Viterbi = np.ones(shape=(len(S)-1, len(text)+1))*-float('inf')\n",
    "    pointers = [[False]+[False for t in text] for s in S[:-1]] \n",
    "    Viterbi[0][0] = np.log(1)\n",
    "    pointers[0][0] = 'root'\n",
    "    offsets = {'M':(3,1),'D':(4,0), 'I':(2,1)}\n",
    "\n",
    "    # init 1st column all del states (not match or ins)\n",
    "    prev = (0)\n",
    "    for state in range(len(S)):\n",
    "        if S[state][0]=='D':\n",
    "            Viterbi[state][0] = Viterbi[prev][0] + np.log(T[prev][state])\n",
    "            pointers[state][0] = (prev,0)\n",
    "            prev=state\n",
    "\n",
    "    for column in range(1, len(text)+1):\n",
    "        symbol = alpha.index(text[column-1])\n",
    "        for state in range(1, len(S[:-1])):\n",
    "            (x,y) = offsets[S[state][0]]\n",
    "            best_score = -float('inf')\n",
    "            if S[state][0]=='D':\n",
    "                p_emit = 1\n",
    "            else:\n",
    "                p_emit = E[state][symbol]\n",
    "#             print('\\nCurrent cell', S[state],(state, column), \"; symbol:\", text[column-1],\n",
    "#                   \"p_emit =\", np.around(p_emit,4))\n",
    "            for prev in range(state-x, state-x+3):\n",
    "                if not prev<0:\n",
    "                    p_prev = Viterbi[prev][column-y]\n",
    "                    p_trans = T[prev][state]\n",
    "                    p_edge = p_prev + np.log(p_trans)\n",
    "#                     print(\"\\tedge:\",S[prev], (prev, column-y), '–>',  S[state], (state, column),)\n",
    "#                     print(\"\\t\\tViterbi: {}, Transition: {};;; Edge -> {}\".format(np.around(p_prev,4),\n",
    "#                                                                                  np.around(p_trans,4),\n",
    "#                                                                                  np.around(p_edge,4)))\n",
    "                    if p_edge > best_score:\n",
    "                        best_score = p_edge\n",
    "                        Viterbi[state][column] = p_edge + np.log(p_emit)\n",
    "                        pointers[state][column] = (prev, column-y)\n",
    "    #                     print(\"\\t\\t\\t**best edge:\", (prev, column-y), S[prev])\n",
    "\n",
    "    print(pd.DataFrame(np.around(Viterbi,3), columns=['start']+list(text), index=S[:-1]))\n",
    "    print(pd.DataFrame(pointers, columns=['start']+list(text), index=S[:-1]))                \n",
    "    return Viterbi, pointers\n",
    "\n",
    "def backtrack(Viterbi, pointers):\n",
    "    ### walk back pointers:\n",
    "    # start backtrack from greatest probability in last column of viterbi\n",
    "    score = -float('inf')\n",
    "    for state in range(len(S)-4,len(S)-1):\n",
    "        if Viterbi[state][len(text)] > score:\n",
    "            last = state\n",
    "            score = Viterbi[state][len(text)]\n",
    "    path = ['']\n",
    "    print(\"Best Score at end of Path:\", score, \"at state:\", S[last])\n",
    "\n",
    "    # backtrack to recreate max likelihood hidden_path in reverse\n",
    "    x,y = last, len(text)\n",
    "    while path:\n",
    "        path.append(S[x])    \n",
    "        (x,y) = pointers[x][y]\n",
    "        if (x,y) ==(0,0):\n",
    "            path = ' '.join(path[-1:0:-1])\n",
    "            break\n",
    "\n",
    "    print('\\n\\nPATH:\\n',path)\n",
    "    return path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text:  AEFDFDC\n",
      "Alignment:\n",
      "\t ACDEFACADF\n",
      "\t AFDA---CCF\n",
      "\t A--EFD-FDC\n",
      "\t ACAEF--A-C\n",
      "\t ADDEFAAADF\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4, 'F': 5}\n",
      "Theta: 0.4\n",
      "Sigma: 0.01\n",
      "      S     I0     M1     D1     I1     M2     D2     I2     M3     D3  ...  \\\n",
      "S   0.0  0.010  0.981  0.010  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I0  0.0  0.333  0.333  0.333  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M1  0.0  0.000  0.000  0.000  0.010  0.786  0.204  0.000  0.000  0.000  ...   \n",
      "D1  0.0  0.000  0.000  0.000  0.333  0.333  0.333  0.000  0.000  0.000  ...   \n",
      "I1  0.0  0.000  0.000  0.000  0.333  0.333  0.333  0.000  0.000  0.000  ...   \n",
      "M2  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.010  0.981  0.010  ...   \n",
      "D2  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.010  0.010  0.981  ...   \n",
      "I2  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.333  0.333  0.333  ...   \n",
      "M3  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D3  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I3  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M4  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D4  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I4  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M5  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D5  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I5  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M6  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D6  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I6  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M7  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D7  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I7  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "M8  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "D8  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "I8  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "E   0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  ...   \n",
      "\n",
      "       M6    D6     I6     M7     D7     I7     M8     D8    I8     E  \n",
      "S   0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I0  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M1  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "D1  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I1  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M2  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "D2  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I2  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M3  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "D3  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I3  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M4  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "D4  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I4  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M5  0.252  0.01  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "D5  0.981  0.01  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "I5  0.592  0.01  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "M6  0.000  0.00  0.010  0.786  0.204  0.000  0.000  0.000  0.00  0.00  \n",
      "D6  0.000  0.00  0.333  0.333  0.333  0.000  0.000  0.000  0.00  0.00  \n",
      "I6  0.000  0.00  0.333  0.333  0.333  0.000  0.000  0.000  0.00  0.00  \n",
      "M7  0.000  0.00  0.000  0.000  0.000  0.010  0.981  0.010  0.00  0.00  \n",
      "D7  0.000  0.00  0.000  0.000  0.000  0.010  0.981  0.010  0.00  0.00  \n",
      "I7  0.000  0.00  0.000  0.000  0.000  0.333  0.333  0.333  0.00  0.00  \n",
      "M8  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.01  0.99  \n",
      "D8  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.50  0.50  \n",
      "I8  0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.50  0.50  \n",
      "E   0.000  0.00  0.000  0.000  0.000  0.000  0.000  0.000  0.00  0.00  \n",
      "\n",
      "[27 rows x 27 columns]\n",
      "        A      B      C      D      E      F\n",
      "S   0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I0  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M1  0.953  0.009  0.009  0.009  0.009  0.009\n",
      "D1  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I1  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M2  0.009  0.009  0.481  0.245  0.009  0.245\n",
      "D2  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I2  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M3  0.245  0.009  0.009  0.717  0.009  0.009\n",
      "D3  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I3  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M4  0.198  0.009  0.009  0.009  0.764  0.009\n",
      "D4  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I4  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M5  0.009  0.009  0.009  0.009  0.009  0.953\n",
      "D5  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I5  0.575  0.009  0.198  0.198  0.009  0.009\n",
      "M6  0.575  0.009  0.198  0.009  0.009  0.198\n",
      "D6  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I6  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M7  0.009  0.009  0.245  0.717  0.009  0.009\n",
      "D7  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I7  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "M8  0.009  0.009  0.387  0.009  0.009  0.575\n",
      "D8  0.000  0.000  0.000  0.000  0.000  0.000\n",
      "I8  0.167  0.167  0.167  0.167  0.167  0.167\n",
      "E   0.000  0.000  0.000  0.000  0.000  0.000\n",
      "     start       A       E       F       D       F       D       C\n",
      "S    0.000    -inf    -inf    -inf    -inf    -inf    -inf    -inf\n",
      "I0    -inf  -6.426  -9.317 -12.207 -15.098 -17.988 -20.878 -23.769\n",
      "M1    -inf  -0.068 -12.189 -15.079 -17.969 -20.860 -23.750 -26.640\n",
      "D1  -4.635  -7.525 -10.415 -13.306 -16.196 -19.087 -21.977 -24.867\n",
      "I1    -inf  -7.525  -6.494  -9.385 -12.275 -15.166 -18.056 -20.946\n",
      "M2    -inf -10.397  -4.972  -8.998 -11.889 -14.779 -17.669 -19.886\n",
      "D2  -5.733  -1.658  -7.593 -10.483 -13.374 -16.264 -19.155 -22.045\n",
      "I2    -inf -12.160  -8.085 -10.975 -13.865 -16.756 -19.646 -22.536\n",
      "M3    -inf -11.773 -10.956  -9.655  -9.351 -16.572 -15.131 -22.353\n",
      "D3  -5.753  -1.678  -7.613 -10.503 -13.393 -16.284 -19.174 -22.064\n",
      "I3    -inf -12.179  -8.104 -10.995 -13.885 -15.777 -18.668 -21.558\n",
      "M4    -inf  -7.391  -1.966 -12.296 -14.338 -14.034 -20.967 -19.814\n",
      "D4 -10.388  -6.312  -9.203 -12.093 -13.985 -16.876 -19.766 -22.657\n",
      "I4    -inf -13.278  -9.203  -8.393 -11.283 -14.174 -17.064 -19.954\n",
      "M5    -inf -16.150 -12.075  -2.255 -14.155 -12.430 -18.937 -22.826\n",
      "D5 -11.486  -7.411  -3.557  -9.491 -12.382 -15.272 -18.163 -21.053\n",
      "I5    -inf -16.674 -16.709 -12.855  -4.178  -9.762 -12.303 -14.843\n",
      "M6    -inf -12.058 -12.094  -5.195  -8.295  -6.321 -14.950 -14.445\n",
      "D6 -16.121 -12.046  -8.191  -6.890  -8.813 -14.397 -16.937 -19.477\n",
      "I6    -inf -19.011 -14.936 -11.082  -9.780 -11.703 -12.747 -15.637\n",
      "M7    -inf -21.883 -16.962 -13.953  -5.768 -13.199  -6.894 -15.251\n",
      "D7 -17.220 -13.144  -9.290  -6.785  -9.885  -7.911 -13.846 -16.035\n",
      "I7    -inf -23.646 -19.571 -15.716 -13.212 -12.195 -14.337 -13.320\n",
      "M8    -inf -21.903 -17.827  -9.862 -11.468  -6.340 -12.594  -7.863\n",
      "D8 -21.854 -17.779 -13.925 -11.420 -10.403 -12.546 -11.528 -14.419\n",
      "I8    -inf -24.339 -20.264 -16.410 -13.905 -12.888 -12.757 -14.013\n",
      "      start        A        E        F        D        F        D        C\n",
      "S      root    False    False    False    False    False    False    False\n",
      "I0    False   (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)\n",
      "M1    False   (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)\n",
      "D1   (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)   (1, 7)\n",
      "I1    False   (3, 0)   (2, 1)   (4, 2)   (4, 3)   (4, 4)   (4, 5)   (4, 6)\n",
      "M2    False   (3, 0)   (2, 1)   (4, 2)   (4, 3)   (4, 4)   (4, 5)   (4, 6)\n",
      "D2   (3, 0)   (2, 1)   (4, 2)   (4, 3)   (4, 4)   (4, 5)   (4, 6)   (4, 7)\n",
      "I2    False   (6, 0)   (6, 1)   (7, 2)   (7, 3)   (7, 4)   (7, 5)   (7, 6)\n",
      "M3    False   (6, 0)   (6, 1)   (5, 2)   (5, 3)   (5, 4)   (5, 5)   (5, 6)\n",
      "D3   (6, 0)   (6, 1)   (6, 2)   (6, 3)   (6, 4)   (6, 5)   (6, 6)   (6, 7)\n",
      "I3    False   (9, 0)   (9, 1)  (10, 2)  (10, 3)   (8, 4)  (10, 5)   (8, 6)\n",
      "M4    False   (9, 0)   (9, 1)   (9, 2)   (8, 3)   (8, 4)   (9, 5)   (8, 6)\n",
      "D4   (9, 0)   (9, 1)  (10, 2)  (10, 3)   (8, 4)  (10, 5)   (8, 6)  (10, 7)\n",
      "I4    False  (12, 0)  (12, 1)  (11, 2)  (13, 3)  (13, 4)  (13, 5)  (13, 6)\n",
      "M5    False  (12, 0)  (12, 1)  (11, 2)  (13, 3)  (13, 4)  (11, 5)  (13, 6)\n",
      "D5  (12, 0)  (12, 1)  (11, 2)  (13, 3)  (13, 4)  (13, 5)  (13, 6)  (13, 7)\n",
      "I5    False  (15, 0)  (15, 1)  (15, 2)  (14, 3)  (16, 4)  (16, 5)  (16, 6)\n",
      "M6    False  (15, 0)  (15, 1)  (15, 2)  (14, 3)  (16, 4)  (16, 5)  (16, 6)\n",
      "D6  (15, 0)  (15, 1)  (15, 2)  (14, 3)  (16, 4)  (16, 5)  (16, 6)  (16, 7)\n",
      "I6    False  (18, 0)  (18, 1)  (18, 2)  (18, 3)  (18, 4)  (17, 5)  (19, 6)\n",
      "M7    False  (18, 0)  (17, 1)  (18, 2)  (17, 3)  (17, 4)  (17, 5)  (19, 6)\n",
      "D7  (18, 0)  (18, 1)  (18, 2)  (17, 3)  (17, 4)  (17, 5)  (19, 6)  (17, 7)\n",
      "I7    False  (21, 0)  (21, 1)  (21, 2)  (21, 3)  (20, 4)  (21, 5)  (20, 6)\n",
      "M8    False  (21, 0)  (21, 1)  (21, 2)  (21, 3)  (20, 4)  (21, 5)  (20, 6)\n",
      "D8  (21, 0)  (21, 1)  (21, 2)  (21, 3)  (20, 4)  (21, 5)  (20, 6)  (22, 7)\n",
      "I8    False  (24, 0)  (24, 1)  (24, 2)  (24, 3)  (24, 4)  (23, 5)  (24, 6)\n",
      "Best Score at end of Path: -7.863084063165768 at state: M8\n",
      "\n",
      "\n",
      "PATH:\n",
      " M1 D2 D3 M4 M5 I5 M6 M7 M8\n"
     ]
    }
   ],
   "source": [
    "### Sample data:\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/ba10g.sample.txt') as f:\n",
    "    text, alignment, alpha, theta, sigma = parse_HMM(f)\n",
    "\n",
    "S, T, E, seed = profileHMM(alignment, alpha, theta)\n",
    "alpha = list(alpha.keys())\n",
    "print(pd.DataFrame(np.around(T,3), index = S, columns=S))\n",
    "print(pd.DataFrame(np.around(E,3), index = S, columns=alpha))\n",
    "V, P = Viterbi_alignment(T,E,S, alpha, text)\n",
    "path = backtrack(V,P)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting the right answer on the sample dataset now! yay! going to call it quits for tonight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text:  EEBEABDCEEABCCCEEBDEDCADEDACCDCBBEECDBDACABDADCBE\n",
      "Alignment:\n",
      "\t EEBBA--C-DBAA-AECD--BDB---CC-DDCBBCEDE-EBB-DAEE-C\n",
      "\t EEEEABBCEABBCDEE-DAEBDBAEDC-BDBCB--C-B-BCA-DAEECA\n",
      "\t --CEB-ACCDEACEEEEDBEED-ADBCCDAC--BC--BDBCAEDAEECC\n",
      "\t A--AABDCE-A-CD-ECD-EBBA-EDC-DACCBBCCD-D-BA-DAAEBC\n",
      "\t EECEAB--EDDACCE-CD-E--B-EDCCD-CCBBCCD-DBBA--AE-CA\n",
      "\t E-CDA-DCECAAECB-CDCEB-B-BDCCD---B--CD-DBCDBDAEB-C\n",
      "\t EBCEAEDC-DABC--A-DCEDDBAED-CD-CCBBCCEBDB--BEA-EEC\n",
      "\t AC-E-BDCEDAADDEECDEEB-CAEDC-DD-CBBCCD-DBCABDAEECC\n",
      "\t EECBABDCEDEAEC-DCDC-BDBDEDA-D-AD-A-EABEB--BDA-ECC\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4}\n",
      "Theta: 0.359\n",
      "Sigma: 0.01\n",
      "       S     I0     M1     D1     I1     M2     D2   I2   M3   D3  ...  M42  \\\n",
      "S    0.0  0.010  0.873  0.118  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I0   0.0  0.333  0.333  0.333  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M1   0.0  0.000  0.000  0.000  0.010  0.738  0.252  0.0  0.0  0.0  ...  0.0   \n",
      "D1   0.0  0.000  0.000  0.000  0.010  0.010  0.981  0.0  0.0  0.0  ...  0.0   \n",
      "I1   0.0  0.000  0.000  0.000  0.333  0.333  0.333  0.0  0.0  0.0  ...  0.0   \n",
      "..   ...    ...    ...    ...    ...    ...    ...  ...  ...  ...  ...  ...   \n",
      "I43  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M44  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "D44  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I44  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "E    0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "\n",
      "     D42  I42  M43  D43    I43    M44    D44   I44     E  \n",
      "S    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I0   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "M1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "D1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "..   ...  ...  ...  ...    ...    ...    ...   ...   ...  \n",
      "I43  0.0  0.0  0.0  0.0  0.333  0.333  0.333  0.00  0.00  \n",
      "M44  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.01  0.99  \n",
      "D44  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.50  0.50  \n",
      "I44  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.50  0.50  \n",
      "E    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "\n",
      "[135 rows x 135 columns]\n",
      "         A     B     C     D      E\n",
      "S    0.000  0.00  0.00  0.00  0.000\n",
      "I0   0.200  0.20  0.20  0.20  0.200\n",
      "M1   0.248  0.01  0.01  0.01  0.724\n",
      "D1   0.000  0.00  0.00  0.00  0.000\n",
      "I1   0.200  0.20  0.20  0.20  0.200\n",
      "..     ...   ...   ...   ...    ...\n",
      "I43  0.200  0.20  0.20  0.20  0.200\n",
      "M44  0.221  0.01  0.75  0.01  0.010\n",
      "D44  0.000  0.00  0.00  0.00  0.000\n",
      "I44  0.200  0.20  0.20  0.20  0.200\n",
      "E    0.000  0.00  0.00  0.00  0.000\n",
      "\n",
      "[135 rows x 5 columns]\n",
      "       start        E        E        B        E        A        B       D  \\\n",
      "S      0.000     -inf     -inf     -inf     -inf     -inf     -inf    -inf   \n",
      "I0      -inf   -6.244   -8.952  -11.660  -14.368  -17.076  -19.784 -22.492   \n",
      "M1      -inf   -0.459   -7.666  -14.705  -13.082  -16.863  -22.829 -25.537   \n",
      "D1    -2.141   -7.343  -10.051  -12.759  -15.467  -18.175  -20.883 -23.591   \n",
      "I1      -inf   -8.385   -6.704   -9.412  -12.120  -14.828  -17.536 -20.244   \n",
      "..       ...      ...      ...      ...      ...      ...      ...     ...   \n",
      "D43 -125.092 -117.390 -109.759 -102.609  -94.906  -93.409  -87.587 -89.204   \n",
      "I43     -inf -131.336 -123.634 -116.004 -108.853 -101.151  -99.653 -93.831   \n",
      "M44     -inf -129.766 -122.064 -114.433 -107.282  -96.435  -98.083 -92.261   \n",
      "D44 -129.727 -122.025 -114.394 -107.243  -99.541  -98.044  -92.222 -93.839   \n",
      "I44     -inf -132.030 -124.327 -116.697 -109.546 -101.844 -100.347 -94.525   \n",
      "\n",
      "          C       E  ...        A        C        A        B        D  \\\n",
      "S      -inf    -inf  ...     -inf     -inf     -inf     -inf     -inf   \n",
      "I0  -25.201 -27.909  ... -111.858 -114.566 -117.274 -119.982 -122.690   \n",
      "M1  -28.245 -26.622  ... -111.645 -117.611 -117.061 -123.027 -125.735   \n",
      "D1  -26.299 -29.007  ... -112.957 -115.665 -118.373 -121.081 -123.789   \n",
      "I1  -22.952 -25.660  ... -109.609 -112.318 -115.026 -117.734 -120.442   \n",
      "..      ...     ...  ...      ...      ...      ...      ...      ...   \n",
      "D43 -89.375 -82.880  ...  -79.298  -76.079  -79.664  -81.409  -86.982   \n",
      "I43 -95.448 -94.595  ...  -81.209  -83.917  -81.298  -84.007  -86.715   \n",
      "M44 -89.511 -93.024  ...  -81.109  -79.605  -76.583  -84.334  -85.192   \n",
      "D44 -92.985 -87.514  ...  -82.308  -79.689  -82.397  -85.105  -87.813   \n",
      "I44 -96.142 -95.288  ...  -80.686  -82.989  -81.992  -82.817  -85.120   \n",
      "\n",
      "           A        D        C        B        E  \n",
      "S       -inf     -inf     -inf     -inf     -inf  \n",
      "I0  -125.398 -128.106 -130.814 -133.523 -136.231  \n",
      "M1  -125.185 -131.151 -133.859 -136.567 -134.944  \n",
      "D1  -126.497 -129.205 -131.913 -134.621 -137.329  \n",
      "I1  -123.150 -125.858 -128.566 -131.274 -133.982  \n",
      "..       ...      ...      ...      ...      ...  \n",
      "D43  -84.982  -81.763  -85.348  -87.721  -88.549  \n",
      "I43  -89.423  -91.226  -88.007  -87.306  -90.014  \n",
      "M44  -86.519  -89.655  -82.070  -85.735  -90.876  \n",
      "D44  -89.617  -86.397  -85.696  -88.404  -91.112  \n",
      "I44  -87.422  -89.725  -88.700  -87.999  -90.301  \n",
      "\n",
      "[134 rows x 50 columns]\n",
      "        start         E         E         B         E         A         B  \\\n",
      "S        root     False     False     False     False     False     False   \n",
      "I0      False    (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)   \n",
      "M1      False    (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)   \n",
      "D1     (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)    (1, 6)   \n",
      "I1      False    (3, 0)    (2, 1)    (4, 2)    (4, 3)    (4, 4)    (4, 5)   \n",
      "..        ...       ...       ...       ...       ...       ...       ...   \n",
      "D43  (126, 0)  (125, 1)  (125, 2)  (126, 3)  (125, 4)  (126, 5)  (125, 6)   \n",
      "I43     False  (129, 0)  (129, 1)  (129, 2)  (129, 3)  (129, 4)  (129, 5)   \n",
      "M44     False  (129, 0)  (129, 1)  (129, 2)  (129, 3)  (129, 4)  (129, 5)   \n",
      "D44  (129, 0)  (129, 1)  (129, 2)  (129, 3)  (129, 4)  (129, 5)  (129, 6)   \n",
      "I44     False  (132, 0)  (132, 1)  (132, 2)  (132, 3)  (132, 4)  (132, 5)   \n",
      "\n",
      "            D         C         E  ...          A          C          A  \\\n",
      "S       False     False     False  ...      False      False      False   \n",
      "I0     (1, 6)    (1, 7)    (1, 8)  ...    (1, 39)    (1, 40)    (1, 41)   \n",
      "M1     (1, 6)    (1, 7)    (1, 8)  ...    (1, 39)    (1, 40)    (1, 41)   \n",
      "D1     (1, 7)    (1, 8)    (1, 9)  ...    (1, 40)    (1, 41)    (1, 42)   \n",
      "I1     (4, 6)    (4, 7)    (4, 8)  ...    (4, 39)    (4, 40)    (4, 41)   \n",
      "..        ...       ...       ...  ...        ...        ...        ...   \n",
      "D43  (125, 7)  (125, 8)  (125, 9)  ...  (126, 40)  (125, 41)  (125, 42)   \n",
      "I43  (129, 6)  (129, 7)  (128, 8)  ...  (130, 39)  (130, 40)  (128, 41)   \n",
      "M44  (129, 6)  (129, 7)  (128, 8)  ...  (130, 39)  (129, 40)  (128, 41)   \n",
      "D44  (129, 7)  (128, 8)  (129, 9)  ...  (130, 40)  (128, 41)  (130, 42)   \n",
      "I44  (132, 6)  (132, 7)  (132, 8)  ...  (133, 39)  (133, 40)  (132, 41)   \n",
      "\n",
      "             B          D          A          D          C          B  \\\n",
      "S        False      False      False      False      False      False   \n",
      "I0     (1, 42)    (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)   \n",
      "M1     (1, 42)    (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)   \n",
      "D1     (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)    (1, 48)   \n",
      "I1     (4, 42)    (4, 43)    (4, 44)    (4, 45)    (4, 46)    (4, 47)   \n",
      "..         ...        ...        ...        ...        ...        ...   \n",
      "D43  (125, 43)  (125, 44)  (126, 45)  (125, 46)  (125, 47)  (125, 48)   \n",
      "I43  (130, 42)  (130, 43)  (130, 44)  (129, 45)  (129, 46)  (128, 47)   \n",
      "M44  (128, 42)  (128, 43)  (128, 44)  (129, 45)  (129, 46)  (128, 47)   \n",
      "D44  (130, 43)  (130, 44)  (129, 45)  (129, 46)  (128, 47)  (130, 48)   \n",
      "I44  (131, 42)  (133, 43)  (133, 44)  (133, 45)  (132, 46)  (132, 47)   \n",
      "\n",
      "             E  \n",
      "S        False  \n",
      "I0     (1, 48)  \n",
      "M1     (1, 48)  \n",
      "D1     (1, 49)  \n",
      "I1     (4, 48)  \n",
      "..         ...  \n",
      "D43  (125, 49)  \n",
      "I43  (130, 48)  \n",
      "M44  (128, 48)  \n",
      "D44  (130, 49)  \n",
      "I44  (133, 48)  \n",
      "\n",
      "[134 rows x 50 columns]\n",
      "Best Score at end of Path: -90.30127149040167 at state: I44\n",
      "\n",
      "\n",
      "PATH:\n",
      " M1 M2 M3 M4 M5 M6 M7 M8 M9 D10 M11 M12 I12 I12 M13 M14 M15 M16 M17 M18 D19 M20 M21 M22 M23 I23 M24 M25 M26 I26 I26 M27 D28 M29 M30 M31 I31 I31 D32 M33 M34 I34 M35 M36 M37 M38 I38 M39 M40 D41 M42 M43 D44 I44 I44\n"
     ]
    }
   ],
   "source": [
    "### extra data:\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/ba10g.extra.txt') as f:\n",
    "    text, alignment, alpha, theta, sigma = parse_HMM(f)\n",
    "\n",
    "\n",
    "S, T, E, seed = profileHMM(alignment, alpha, theta)\n",
    "alpha = list(alpha.keys())\n",
    "print(pd.DataFrame(np.around(T,3), index = S, columns=S))\n",
    "print(pd.DataFrame(np.around(E,3), index = S, columns=alpha))\n",
    "V, P = Viterbi_alignment(T,E,S,alpha, text)\n",
    "path = backtrack(V,P)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TRYING TO FIGURE WHY IT DOESNT WORK ON THIS DATASET????\n",
    "\n",
    "# path = path.split(' ')\n",
    "# len(path)\n",
    "# sol = 'M1 M2 M3 M4 M5 M6 M7 M8 M9 D10 M11 M12 I12 I12 M13 M14 M15 M16 M17 M18 D19 M20 M21 M22 M23 I23 M24 M25 M26 I26 I26 M27 D28 M29 M30 M31 I31 I31 D32 M33 M34 I34 M35 M36 M37 M38 I38 M39 M40 M41 M42 M43 M44'.split(' ')\n",
    "# len(sol)\n",
    "# sum(1 if s[0] in ('M','I') else 0 for s in path), sum(1 if s[0] in ('D') else 0 for s in path)\n",
    "# sum(1 if s[0] in ('M','I') else 0 for s in sol), sum(1 if s[0] in ('D') else 0 for s in sol)\n",
    "# v = pd.DataFrame(np.around(V,3), columns=['start']+list(text), index=S[:-1])\n",
    "# p = pd.DataFrame(P, columns=['start']+list(text), index=S[:-1])     \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text:  ADDBACBDDDBACEABEEBDBEEDAECEAEBDCEADEEBCEADBDCABD\n",
      "Alignment:\n",
      "\t ADD-AEBCDEBECE-BE-BA-EEEAB-ECEBDCEABEE-B--DBEEB-D\n",
      "\t A--EAEB-DEBEC-DB-EBCBEAAAB-EADBDC-AD-EDA-A-B-ABB-\n",
      "\t -DAD-EB-DEBECBDBD-B--EDEABBE-ABD-EADEEBC-ADEDCBBD\n",
      "\t CADCAE-D-E-E--DBEEBA-EEEC-BEADBCCEDABCBC-ADBD--BD\n",
      "\t ADDDA-BED-BE-E-BEEB-BEEDA--EBDB-CEA-ECDCCA-BDABBB\n",
      "\t ADDDAEDDD-BD-EDBEE-DD-E-D-B-ADBD-EB--EBCCE-BDABBD\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4}\n",
      "Theta: 0.327\n",
      "Sigma: 0.01\n",
      "       S     I0     M1     D1     I1     M2     D2   I2   M3   D3  ...  M32  \\\n",
      "S    0.0  0.010  0.819  0.172  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I0   0.0  0.333  0.333  0.333  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M1   0.0  0.000  0.000  0.000  0.010  0.786  0.204  0.0  0.0  0.0  ...  0.0   \n",
      "D1   0.0  0.000  0.000  0.000  0.010  0.981  0.010  0.0  0.0  0.0  ...  0.0   \n",
      "I1   0.0  0.000  0.000  0.000  0.333  0.333  0.333  0.0  0.0  0.0  ...  0.0   \n",
      "..   ...    ...    ...    ...    ...    ...    ...  ...  ...  ...  ...  ...   \n",
      "I33  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M34  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "D34  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I34  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "E    0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "\n",
      "     D32  I32  M33  D33    I33    M34    D34   I34     E  \n",
      "S    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I0   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "M1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "D1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "..   ...  ...  ...  ...    ...    ...    ...   ...   ...  \n",
      "I33  0.0  0.0  0.0  0.0  0.333  0.333  0.333  0.00  0.00  \n",
      "M34  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.01  0.99  \n",
      "D34  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.01  0.99  \n",
      "I34  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.50  0.50  \n",
      "E    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "\n",
      "[105 rows x 105 columns]\n",
      "         A     B     C      D     E\n",
      "S    0.000  0.00  0.00  0.000  0.00\n",
      "I0   0.200  0.20  0.20  0.200  0.20\n",
      "M1   0.771  0.01  0.20  0.010  0.01\n",
      "D1   0.000  0.00  0.00  0.000  0.00\n",
      "I1   0.200  0.20  0.20  0.200  0.20\n",
      "..     ...   ...   ...    ...   ...\n",
      "I33  0.200  0.20  0.20  0.200  0.20\n",
      "M34  0.010  0.20  0.01  0.771  0.01\n",
      "D34  0.000  0.00  0.00  0.000  0.00\n",
      "I34  0.200  0.20  0.20  0.200  0.20\n",
      "E    0.000  0.00  0.00  0.000  0.00\n",
      "\n",
      "[105 rows x 5 columns]\n",
      "       start        A        D        D        B       A       C       B  \\\n",
      "S      0.000     -inf     -inf     -inf     -inf    -inf    -inf    -inf   \n",
      "I0      -inf   -6.244   -8.952  -11.660  -14.368 -17.076 -19.784 -22.492   \n",
      "M1      -inf   -0.459  -11.997  -14.705  -17.413 -15.726 -19.784 -25.537   \n",
      "D1    -1.763   -7.343  -10.051  -12.759  -15.467 -18.175 -20.883 -23.591   \n",
      "I1      -inf   -8.007   -6.704   -9.412  -12.120 -14.828 -17.536 -20.244   \n",
      "..       ...      ...      ...      ...      ...     ...     ...     ...   \n",
      "D33 -105.483  -97.862  -94.018  -91.013  -87.507 -83.435 -79.924 -77.424   \n",
      "I33     -inf -111.727 -104.106 -100.262  -97.257 -92.681 -89.679 -86.168   \n",
      "M34     -inf -110.157  -98.141  -94.297  -92.642 -91.331 -88.109 -81.553   \n",
      "D34 -110.118 -102.497  -98.653  -95.647  -88.027 -88.070 -84.558 -76.938   \n",
      "I34     -inf -116.352 -108.731 -104.376 -100.532 -94.261 -94.304 -90.793   \n",
      "\n",
      "          D       D  ...        C        E        A        D        B  \\\n",
      "S      -inf    -inf  ...     -inf     -inf     -inf     -inf     -inf   \n",
      "I0  -25.201 -27.909  ... -111.858 -114.566 -117.274 -119.982 -122.690   \n",
      "M1  -28.245 -30.953  ... -111.858 -117.611 -115.924 -123.027 -125.735   \n",
      "D1  -26.299 -29.007  ... -112.957 -115.665 -118.373 -121.081 -123.789   \n",
      "I1  -22.952 -25.660  ... -109.610 -112.318 -115.026 -117.734 -120.442   \n",
      "..      ...     ...  ...      ...      ...      ...      ...      ...   \n",
      "D33 -75.282 -75.340  ...  -80.435  -79.000  -79.038  -77.843  -79.225   \n",
      "I33 -81.592 -81.526  ...  -79.401  -82.109  -84.817  -85.283  -84.088   \n",
      "M34 -75.847 -75.561  ...  -78.051  -82.918  -83.673  -79.318  -79.472   \n",
      "D34 -79.696 -76.911  ...  -79.614  -82.064  -80.629  -82.478  -74.857   \n",
      "I34 -83.172 -82.082  ...  -80.981  -83.284  -85.586  -86.863  -85.552   \n",
      "\n",
      "           D        C        A        B        D  \n",
      "S       -inf     -inf     -inf     -inf     -inf  \n",
      "I0  -125.398 -128.106 -130.814 -133.523 -136.231  \n",
      "M1  -128.443 -128.106 -129.465 -136.567 -139.275  \n",
      "D1  -126.497 -129.205 -131.913 -134.621 -137.329  \n",
      "I1  -123.150 -125.858 -128.566 -131.274 -133.982  \n",
      "..       ...      ...      ...      ...      ...  \n",
      "D33  -77.584  -80.292  -79.453  -80.874  -86.626  \n",
      "I33  -79.511  -82.219  -83.867  -85.697  -84.386  \n",
      "M34  -73.767  -82.257  -82.517  -81.082  -78.641  \n",
      "D34  -80.610  -79.213  -84.087  -79.732  -85.484  \n",
      "I34  -81.092  -80.001  -82.304  -84.607  -85.966  \n",
      "\n",
      "[104 rows x 50 columns]\n",
      "       start         A         D         D         B         A         C  \\\n",
      "S       root     False     False     False     False     False     False   \n",
      "I0     False    (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)   \n",
      "M1     False    (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)   \n",
      "D1    (0, 0)    (1, 1)    (1, 2)    (1, 3)    (1, 4)    (1, 5)    (1, 6)   \n",
      "I1     False    (3, 0)    (2, 1)    (4, 2)    (4, 3)    (4, 4)    (4, 5)   \n",
      "..       ...       ...       ...       ...       ...       ...       ...   \n",
      "D33  (96, 0)   (96, 1)   (96, 2)   (96, 3)   (96, 4)   (96, 5)   (96, 6)   \n",
      "I33    False   (99, 0)   (99, 1)   (99, 2)   (99, 3)   (98, 4)   (99, 5)   \n",
      "M34    False   (99, 0)   (99, 1)   (99, 2)   (99, 3)   (98, 4)   (99, 5)   \n",
      "D34  (99, 0)   (99, 1)   (99, 2)   (98, 3)   (98, 4)   (99, 5)   (99, 6)   \n",
      "I34    False  (102, 0)  (102, 1)  (101, 2)  (101, 3)  (102, 4)  (102, 5)   \n",
      "\n",
      "            B         D         D  ...          C          E          A  \\\n",
      "S       False     False     False  ...      False      False      False   \n",
      "I0     (1, 6)    (1, 7)    (1, 8)  ...    (1, 39)    (1, 40)    (1, 41)   \n",
      "M1     (1, 6)    (1, 7)    (1, 8)  ...    (1, 39)    (1, 40)    (1, 41)   \n",
      "D1     (1, 7)    (1, 8)    (1, 9)  ...    (1, 40)    (1, 41)    (1, 42)   \n",
      "I1     (4, 6)    (4, 7)    (4, 8)  ...    (4, 39)    (4, 40)    (4, 41)   \n",
      "..        ...       ...       ...  ...        ...        ...        ...   \n",
      "D33   (95, 7)   (96, 8)   (96, 9)  ...   (96, 40)   (96, 41)   (95, 42)   \n",
      "I33   (99, 6)   (98, 7)   (99, 8)  ...   (98, 39)  (100, 40)  (100, 41)   \n",
      "M34   (99, 6)   (98, 7)   (99, 8)  ...   (98, 39)   (98, 40)   (99, 41)   \n",
      "D34   (98, 7)   (98, 8)   (98, 9)  ...   (98, 40)   (98, 41)   (98, 42)   \n",
      "I34  (102, 6)  (102, 7)  (101, 8)  ...  (102, 39)  (103, 40)  (103, 41)   \n",
      "\n",
      "             D          B          D          C          A          B  \\\n",
      "S        False      False      False      False      False      False   \n",
      "I0     (1, 42)    (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)   \n",
      "M1     (1, 42)    (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)   \n",
      "D1     (1, 43)    (1, 44)    (1, 45)    (1, 46)    (1, 47)    (1, 48)   \n",
      "I1     (4, 42)    (4, 43)    (4, 44)    (4, 45)    (4, 46)    (4, 47)   \n",
      "..         ...        ...        ...        ...        ...        ...   \n",
      "D33   (96, 43)   (95, 44)   (96, 45)   (97, 46)   (95, 47)   (95, 48)   \n",
      "I33   (98, 42)   (99, 43)   (98, 44)  (100, 45)   (98, 46)   (99, 47)   \n",
      "M34   (99, 42)   (99, 43)   (98, 44)   (99, 45)   (98, 46)   (99, 47)   \n",
      "D34   (99, 43)   (98, 44)  (100, 45)   (98, 46)   (98, 47)   (98, 48)   \n",
      "I34  (102, 42)  (101, 43)  (102, 44)  (101, 45)  (103, 46)  (103, 47)   \n",
      "\n",
      "             D  \n",
      "S        False  \n",
      "I0     (1, 48)  \n",
      "M1     (1, 48)  \n",
      "D1     (1, 49)  \n",
      "I1     (4, 48)  \n",
      "..         ...  \n",
      "D33   (97, 49)  \n",
      "I33   (98, 48)  \n",
      "M34   (98, 48)  \n",
      "D34  (100, 49)  \n",
      "I34  (102, 48)  \n",
      "\n",
      "[104 rows x 50 columns]\n",
      "Best Score at end of Path: -78.64139298306095 at state: M34\n",
      "\n",
      "\n",
      "PATH:\n",
      " M1 M2 M3 M4 M5 M6 M7 I7 M8 I8 M9 M10 I10 I10 I10 M11 M12 I12 M13 I13 I13 M14 M15 M16 M17 M18 I18 I18 M19 M20 M21 M22 I22 M23 M24 I24 I24 M25 M26 M27 I27 M28 I28 M29 M30 M31 M32 M33 M34\n"
     ]
    }
   ],
   "source": [
    "### stepik test data:\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/dataset_26259_14.txt') as f:\n",
    "    text, alignment, alpha, theta, sigma = parse_HMM(f)\n",
    "\n",
    "\n",
    "S, T, E, seed = profileHMM(alignment, alpha, theta)\n",
    "alpha = list(alpha.keys())\n",
    "print(pd.DataFrame(np.around(T,3), index = S, columns=S))\n",
    "print(pd.DataFrame(np.around(E,3), index = S, columns=alpha))\n",
    "V, P = Viterbi_alignment(T,E,S,alpha,text)\n",
    "path = backtrack(V,P)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### correct on the stepik dataset, now try rosalind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Text:  EACCEECBECCADDEBABEEBEABDCDEEBECEBECBBDCDCDDCDCAA\n",
      "Alignment:\n",
      "\t EACCEEB-A-EAABEBA--ECCDBD-DCACE--BACAB-CB-DABA-D-\n",
      "\t EAC--EBCE--A--EB-DEBCCDB--CCACEBE--DABCADCCABA-AC\n",
      "\t EACCEEBBECEAABEBCEECCC-C-DDC-EECEBBAABACDC-ABAC-C\n",
      "\t EACC-BB-E--A--EBA-EE--D-DCDCCCDB-BD-ABECDCD-BAEA-\n",
      "\t -ACC-E-BEC---B-BA---CED-DCDCBCECEBC-ABACD-D-BA--C\n",
      "\t CACEEE---CEA-BE--DEECBDBD-DCA--B----DB-CDCEBBA-AC\n",
      "\t EACD--BDEC-AAB-BADEE-CD-BCDCAC-DABDCABA-ECDC-AEAE\n",
      "\t EACC--BBECEAABDBADEECCDB-CADACB-E--CAB-C--DAB-BAC\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4}\n",
      "Theta: 0.297\n",
      "Sigma: 0.01\n",
      "       S     I0     M1     D1     I1     M2     D2   I2   M3   D3  ...  M31  \\\n",
      "S    0.0  0.010  0.859  0.131  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I0   0.0  0.333  0.333  0.333  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M1   0.0  0.000  0.000  0.000  0.010  0.981  0.010  0.0  0.0  0.0  ...  0.0   \n",
      "D1   0.0  0.000  0.000  0.000  0.010  0.981  0.010  0.0  0.0  0.0  ...  0.0   \n",
      "I1   0.0  0.000  0.000  0.000  0.333  0.333  0.333  0.0  0.0  0.0  ...  0.0   \n",
      "..   ...    ...    ...    ...    ...    ...    ...  ...  ...  ...  ...  ...   \n",
      "I32  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "M33  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "D33  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "I33  0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "E    0.0  0.000  0.000  0.000  0.000  0.000  0.000  0.0  0.0  0.0  ...  0.0   \n",
      "\n",
      "     D31  I31  M32  D32    I32    M33    D33   I33     E  \n",
      "S    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I0   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "M1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "D1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "I1   0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "..   ...  ...  ...  ...    ...    ...    ...   ...   ...  \n",
      "I32  0.0  0.0  0.0  0.0  0.333  0.333  0.333  0.00  0.00  \n",
      "M33  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.01  0.99  \n",
      "D33  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.01  0.99  \n",
      "I33  0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.50  0.50  \n",
      "E    0.0  0.0  0.0  0.0  0.000  0.000  0.000  0.00  0.00  \n",
      "\n",
      "[102 rows x 102 columns]\n",
      "        A     B      C     D      E\n",
      "S    0.00  0.00  0.000  0.00  0.000\n",
      "I0   0.20  0.20  0.200  0.20  0.200\n",
      "M1   0.01  0.01  0.146  0.01  0.826\n",
      "D1   0.00  0.00  0.000  0.00  0.000\n",
      "I1   0.20  0.20  0.200  0.20  0.200\n",
      "..    ...   ...    ...   ...    ...\n",
      "I32  0.20  0.20  0.200  0.20  0.200\n",
      "M33  0.01  0.01  0.803  0.01  0.168\n",
      "D33  0.00  0.00  0.000  0.00  0.000\n",
      "I33  0.20  0.20  0.200  0.20  0.200\n",
      "E    0.00  0.00  0.000  0.00  0.000\n",
      "\n",
      "[102 rows x 5 columns]\n",
      "       start        E        A       C       C       E       E       C  \\\n",
      "S      0.000     -inf     -inf    -inf    -inf    -inf    -inf    -inf   \n",
      "I0      -inf   -6.244   -8.952 -11.660 -14.368 -17.076 -19.784 -22.492   \n",
      "M1      -inf   -0.343  -11.997 -11.978 -14.686 -15.658 -18.366 -22.810   \n",
      "D1    -2.032   -7.343  -10.051 -12.759 -15.467 -18.175 -20.883 -23.591   \n",
      "I1      -inf   -8.276   -6.587  -9.295 -12.003 -14.711 -17.419 -20.127   \n",
      "..       ...      ...      ...     ...     ...     ...     ...     ...   \n",
      "D32 -103.466  -95.534  -88.231 -85.799 -82.465 -78.347 -75.659 -72.715   \n",
      "I32     -inf -109.710 -101.778 -94.475 -92.043 -88.709 -84.591 -81.903   \n",
      "M33     -inf -105.268 -100.208 -88.470 -86.037 -84.267 -80.149 -75.898   \n",
      "D33 -108.101 -100.169  -92.866 -90.433 -87.099 -82.982 -80.294 -77.349   \n",
      "I33     -inf -114.335 -106.403 -99.100 -94.704 -92.272 -89.216 -86.383   \n",
      "\n",
      "          B       E  ...        C        D        C        D        D  \\\n",
      "S      -inf    -inf  ...     -inf     -inf     -inf     -inf     -inf   \n",
      "I0  -25.201 -27.909  ... -111.858 -114.566 -117.274 -119.982 -122.690   \n",
      "M1  -28.245 -26.490  ... -112.176 -117.611 -117.592 -123.027 -125.735   \n",
      "D1  -26.299 -29.007  ... -112.957 -115.665 -118.373 -121.081 -123.789   \n",
      "I1  -22.836 -25.544  ... -109.493 -112.201 -114.909 -117.617 -120.325   \n",
      "..      ...     ...  ...      ...      ...      ...      ...      ...   \n",
      "D32 -65.411 -62.895  ...  -79.885  -81.741  -82.950  -84.457  -87.700   \n",
      "I32 -78.959 -71.656  ...  -82.945  -85.653  -86.839  -89.194  -89.904   \n",
      "M33 -77.388 -67.213  ...  -77.340  -84.559  -81.234  -87.624  -88.734   \n",
      "D33 -70.046 -67.530  ...  -82.865  -81.693  -86.439  -84.758  -86.284   \n",
      "I33 -82.132 -76.280  ...  -81.340  -83.575  -85.877  -87.468  -89.771   \n",
      "\n",
      "           C        D        C        A        A  \n",
      "S       -inf     -inf     -inf     -inf     -inf  \n",
      "I0  -125.398 -128.106 -130.814 -133.523 -136.231  \n",
      "M1  -125.716 -131.151 -131.132 -136.567 -139.275  \n",
      "D1  -126.497 -129.205 -131.913 -134.621 -137.329  \n",
      "I1  -123.033 -125.741 -128.449 -131.158 -133.866  \n",
      "..       ...      ...      ...      ...      ...  \n",
      "D32  -86.585  -92.402  -93.611  -91.270  -90.841  \n",
      "I32  -91.430  -92.830  -93.539  -96.247  -97.514  \n",
      "M33  -85.825  -91.259  -87.934  -98.285  -95.944  \n",
      "D33  -91.220  -88.394  -94.638  -93.857  -91.534  \n",
      "I33  -92.073  -92.059  -94.362  -94.169  -96.471  \n",
      "\n",
      "[101 rows x 50 columns]\n",
      "       start        E        A        C        C        E        E        C  \\\n",
      "S       root    False    False    False    False    False    False    False   \n",
      "I0     False   (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)   \n",
      "M1     False   (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)   \n",
      "D1    (0, 0)   (1, 1)   (1, 2)   (1, 3)   (1, 4)   (1, 5)   (1, 6)   (1, 7)   \n",
      "I1     False   (3, 0)   (2, 1)   (4, 2)   (4, 3)   (4, 4)   (4, 5)   (4, 6)   \n",
      "..       ...      ...      ...      ...      ...      ...      ...      ...   \n",
      "D32  (93, 0)  (93, 1)  (92, 2)  (92, 3)  (93, 4)  (93, 5)  (92, 6)  (93, 7)   \n",
      "I32    False  (96, 0)  (96, 1)  (96, 2)  (96, 3)  (96, 4)  (96, 5)  (96, 6)   \n",
      "M33    False  (96, 0)  (96, 1)  (96, 2)  (96, 3)  (96, 4)  (96, 5)  (96, 6)   \n",
      "D33  (96, 0)  (96, 1)  (96, 2)  (96, 3)  (96, 4)  (96, 5)  (96, 6)  (96, 7)   \n",
      "I33    False  (99, 0)  (99, 1)  (99, 2)  (98, 3)  (98, 4)  (99, 5)  (98, 6)   \n",
      "\n",
      "           B        E  ...          C         D          C         D  \\\n",
      "S      False    False  ...      False     False      False     False   \n",
      "I0    (1, 7)   (1, 8)  ...    (1, 39)   (1, 40)    (1, 41)   (1, 42)   \n",
      "M1    (1, 7)   (1, 8)  ...    (1, 39)   (1, 40)    (1, 41)   (1, 42)   \n",
      "D1    (1, 8)   (1, 9)  ...    (1, 40)   (1, 41)    (1, 42)   (1, 43)   \n",
      "I1    (4, 7)   (4, 8)  ...    (4, 39)   (4, 40)    (4, 41)   (4, 42)   \n",
      "..       ...      ...  ...        ...       ...        ...       ...   \n",
      "D32  (93, 8)  (94, 9)  ...   (94, 40)  (92, 41)   (94, 42)  (92, 43)   \n",
      "I32  (96, 7)  (96, 8)  ...   (95, 39)  (97, 40)   (95, 41)  (96, 42)   \n",
      "M33  (96, 7)  (96, 8)  ...   (95, 39)  (96, 40)   (95, 41)  (96, 42)   \n",
      "D33  (96, 8)  (96, 9)  ...   (95, 40)  (95, 41)   (95, 42)  (95, 43)   \n",
      "I33  (98, 7)  (99, 8)  ...  (100, 39)  (98, 40)  (100, 41)  (98, 42)   \n",
      "\n",
      "             D          C         D          C         A          A  \n",
      "S        False      False     False      False     False      False  \n",
      "I0     (1, 43)    (1, 44)   (1, 45)    (1, 46)   (1, 47)    (1, 48)  \n",
      "M1     (1, 43)    (1, 44)   (1, 45)    (1, 46)   (1, 47)    (1, 48)  \n",
      "D1     (1, 44)    (1, 45)   (1, 46)    (1, 47)   (1, 48)    (1, 49)  \n",
      "I1     (4, 43)    (4, 44)   (4, 45)    (4, 46)   (4, 47)    (4, 48)  \n",
      "..         ...        ...       ...        ...       ...        ...  \n",
      "D32   (92, 44)   (94, 45)  (92, 46)   (94, 47)  (92, 48)   (92, 49)  \n",
      "I32   (95, 43)   (95, 44)  (96, 45)   (95, 46)  (97, 47)   (96, 48)  \n",
      "M33   (95, 43)   (95, 44)  (96, 45)   (95, 46)  (96, 47)   (96, 48)  \n",
      "D33   (95, 44)   (96, 45)  (95, 46)   (97, 47)  (95, 48)   (95, 49)  \n",
      "I33  (100, 43)  (100, 44)  (98, 45)  (100, 46)  (98, 47)  (100, 48)  \n",
      "\n",
      "[101 rows x 50 columns]\n",
      "Best Score at end of Path: -91.53395823364092 at state: D33\n",
      "\n",
      "\n",
      "PATH:\n",
      " M1 M2 M3 M4 I4 M5 M6 I6 M7 I7 I7 M8 I8 M9 M10 M11 M12 I12 M13 M14 M15 M16 M17 I17 M18 M19 M20 M21 M22 M23 I23 I23 I23 I23 I23 I23 I23 I23 I23 I23 I23 I23 M24 D25 D26 M27 I27 M28 M29 D30 M31 M32 D33\n"
     ]
    }
   ],
   "source": [
    "### rosalind test data:\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/rosalind_ba10g.txt') as f:\n",
    "    text, alignment, alpha, theta, sigma = parse_HMM(f)\n",
    "\n",
    "\n",
    "S, T, E, seed = profileHMM(alignment, alpha, theta)\n",
    "alpha = list(alpha.keys())\n",
    "print(pd.DataFrame(np.around(T,3), index = S, columns=S))\n",
    "print(pd.DataFrame(np.around(E,3), index = S, columns=alpha))\n",
    "V, P = Viterbi_alignment(T,E,S,alpha, text)\n",
    "path = backtrack(V,P)\n"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}